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SUMMARY 


The  primary  objective  of  this  project  is  to  develop  arxl  apply  improved  statistical  methodologies  for 
relating  seismic  magnitudes  to  explosion  yields,  treating  both  magnitudes  and  yields  as  uncertain  vari¬ 
ables  and  using  the  censored  yields  in  the  yield  estimation  procedure.  During  the  past  two  years,  our 
major  efforts  have  been 

[A]  Measure  the  teleseismic  P-wave  magnitudes  of  historical  Soviet  explosions  as  well  as  explosions 
from  other  foreign  test  sites  recorded  at  the  optimal  distance  range  from  20°  to  95°. 

[B]  Perform  various  statistical  analyses  of  the  raw  m/,  and  obtain  the  optimal  network  values. 
Conduct  the  maximum-likelihood  magnitude-yield  regressbns  and  analyze  the  source-depth  scal¬ 
ing  relationship. 

[C]  Conduct  a  theoretical  study  to  investigate  relevant  issues.  Improve  and  document  the  statistical  as 
well  as  the  fonward-modeling  tools  currently  in  use. 

Five  technical  reports  were  submitted  during  the  contract  period  (Aug  1991  -  Nov  1993): 

(1)  TGAL-92-05,  “Path-corrected  body-wave  magnitudes  and  yield  estimates  of  Semipalatinsk  explo¬ 
sions". 

(2)  TGAL-92-1 1 .  “Simultaneous  inversion  of  explosion  size  and  path  attenuation  parameter  with  cru¬ 
stal  phases". 

(3)  TGAL-93-06,  "User's  manual  of  FD2:  a  software  package  for  modeling  seismological  problems 
with  2-d/mens/ona/ //near  finite-difference  method”. 

(4)  TGAL-93-07,  “Statistical  characterization  of  rugged  propagation  paths  with  application  to  Rg 
scattering  study". 

(5)  TGAL-93-05,  “Statistical  study  of  Soviet  nuclear  explosions:  lata,  results,  and  software  tools". 

This  final  report,  TGAL-93-05,  summarizes  our  updated  results  obtained  under  Task  [B]  using  the 
data  collected  under  Task  [A],  We  also  give  detailed  descriptions  of  several  key  algorithms  of  our 
software  tools.  Sample  scripts  and  examples  are  furnished  for  these  routines.  The  forward-modeling 
package,  "fd2",  updated  under  Task  [C]  is  documented  in  two  accompanying  reports,  TGAL-93-06  and 
TGAL-93-07. 

Our  database  of  station  mt,  values  based  on  short-period  vertical-component  (SPZ)  recordings  of 
nuclear  explosions  has  been  expanded  to  252  events  located  at  a  variety  of  test  sites.  16,716  carefully 
measured  station  magnitudes,  along  with  10,055  noise  measurements  and  2,004  clipped  measurements, 
were  fed  into  a  maximum-likelihood  inversion  scheme  which  simultaneously  determines  the  event  size 
and  the  station  correction,  as  well  as  the  specific  path  correction  for  each  source-station  pair.  The 
simultaneously-inferred  path  and  station  corrections  are  related  to  kix)wn  geological/geophysical 
features.  Applying  these  path  and  station  corrections  to  the  raw  station  magnitudes  of  any  individual 
explosion  yields  a  systematic  reduction  in  the  fluctuational  variation  of  station  magnitudes  across  the 
whole  network  with  a  reduction  factor  ranging  from  1 .2  to  3  for  all  Soviet  events  in  our  data  set.  Most 
Novaya  Zemlya  events  exhibit  a  variation  reduction  factor  of  2.  With  these  path-correctedfstation- 
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corrected  mi^Pmax).  the  /”b(Pmax)-^b(^fl)INORSARl  bias  between  the  southwest  and  northeast  subre¬ 
gions  of  the  Soviet's  Balapan  test  site  is  assessed  as  0.07  magnitude  unit  [m.u.],  which  is  significantly 
smaller  than  that  of  previous  studies.  This  bias  can  be  further  reduced  somewhat  when  the  based 
on  the  first  motion,  ,  is  used.  First  motion  of  the  initial  short-period  P  waves  also  appears  to  be  a 
very  favorable  source  measure  for  explosions  fired  in  hard  rock  sites  underlain  by  a  stable  mantle  such 
as  Semipalatinsk.  For  example,  based  on  m^P^)  alone  and  without  any  extra  cratering-to-contained 
correction,  the  Balapan  explosion  of  Jan  15,  1965,  is  estimated  to  have  a  yield  of  120  kt.  The  mt^PJ- 
based  yield  estimate  for  the  JVE  event  of  Sep  14, 1988,  is  112  kt.  Between  100  and  150  kt,  the  m^  bias 
between  Eastern  Kazakh  and  NTS  using  our  m^P^ax)  values  is  0.35  m.u.  Along  with  other  software 
tools  developed  under  this  project,  the  explosion  mi,  dataset  is  being  installed  at  CSS. 
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1.  INTRODUCTION 


The  primary  objective  of  this  proje<^  is  to  develop  and  apply  improved  statistical  methodologies  for 
relating  seismic  magnitudes  to  explosion  yields,  treating  both  the  magnitudes  and  yields  as  uncertain 
variables  and  using  censored  yields  in  the  yield  estimation  procedure.  During  the  past  two  years,  our 
major  efforts  have  been 

[A]  Measure  the  teleseismic  P-wave  magnitudes  of  historical  Soviet  explosions  as  well  as  explosions 
from  other  foreign  test  sites  recorded  at  the  optimal  distance  range  from  20*^  to  95’’. 

[B]  Perform  various  statistical  analyses  of  the  raw  m/,  and  obtain  the  optimal  network  rr^  values.  Con¬ 
duct  the  maximum-likelihood  magnitude-yield  regressions  aixl  analyze  the  source-depth  scaling  relation¬ 
ship. 

[C]  Conduct  a  theoretical  study  to  investigate  relevant  issues.  Improve  and  document  the  statistical  as 
well  as  the  fonvard-modeling  tools  currently  in  use. 

This  final  report  summarizes  our  updated  results  obtained  under  Task  [B]  using  the  data  collected 
under  Task  [A].  We  also  present  detailed  descriptions  of  several  key  algorithms  of  our  software  tools 
developed  under  Task  [C].'  Sample  scripts  and  examples  are  furnished  for  these  routines.  That  is,  we 
not  only  present  our  interpretation  of  these  data,  we  also  explain  how  the  analyses  were  carried  out. 
Thus  this  report  is  actually  a  combination  of  a  technical  summary,  a  programmer’s  guide,  and  a  user's 
manual.  Key  routines  covered  in  this  report  are 

[1]  getmag:  a  routine  for  computing  station  magnitudes, 

[2]  emils  (domie):  a  single-event  maximum-likelihood  estimator, 

[3]  mlglm:\he  maximum-likelihood  general  linear  model. 

[4]  geomap:  a  map-plotting  routine, 

[5]  dwisq  (doIsqS):  magnitude-yield  regression  with  uncertain  x  and  y, 

[6]  guessQ:  time-domain  determination  of  Lg  attenuation  coefficient, 

[7]  domle2:  linear  regression  with  censored  y. 

Under  Task  [A]  we  have  accumulated  28,775  carefully-measured  explosion  mt,  values  for  nuclear 
tests  from  a  variety  of  regions,  with  new  data  primarily  from  WWSSN  [World  Wide  Standard  Seismo¬ 
graph  Network]  recordings  of  Soviet  nuclear  tests.  During  the  past  three  years,  our  database  of  station 
mt  values  based  on  short-period  vertical-component  (SPZ)  recordings  of  body  waves  has  been 
expanded  from  112  events  to  252  events  from  a  variety  of  regions  {cf.  Table  1).  It  consists  of  744 
usable  “a"  (/.e. ,  zero-crossing  to  first  peak),  “b"  (/.e.,  first  peak  to  first  trough),  and  “max"  {/.©.,  max 
peak-to-trough  or  trough-to-peak  in  the  first  5  seconds)  event  phases.^  Between  the  distance  range  of 
20°  and  95°,  there  are  16,716  carefully  measured  signals  abng  with  10,055  noisy  measurements  and 
2,004  clipped  measurements.  The  WWSSN  network  is  still  very  valuable,  because  it  provides  data  with 
a  uniform  instrument  response  recorded  over  a  long  time  span  and  with  good  distribution  around  all  test 
sites. 


'  Our  forward-modeling  package,  ftC.  developed  under  Task  [CJ  is  documented  in  an  accompanying  report,  TGAL-93-6. 

*  1 1  "a"  and  1  "b"  phases  are  not  available  (cf.  Table  3),  and  hence  only  744  =  3  x  252  -  12  phases  are  used  in  this 

study. 
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Table  1.  Explosion  ntt,  Database 

01  Jan  90 

31  Dec  92 

Nuclear  Test  Site 

19 

38 

Nevada  Test  Site,  U.S.A. 

6 

6 

Outside  Nevada  Test  Site,  U  S  A. 

3 

3 

Amchitka  Island,  Aleutians,  U.S.A. 

11 

11 

Azgir,  U.S.S.R. 

0 

8 

Orenburg,  U.S.S.R. 

1 

2 

“PNE”,  U.S.S.R. 

0 

14 

Murzhik  (Konystan),  E.  Kazakh 

9 

21 

Degelen  Mountain,  E.  Kazakh 

12 

79 

Balapan  (Shagan  River),  E.  Kazakh 

18 

30 

Northern  Novaya  Zemlya 

6 

6 

Southern  Novaya  Zemlya 

9 

9 

Ahaggar,  French  Sahara 

11 

11 

Tuamoto  Islands,  France 

1 

1 

Rajasthan,  India 

6 

13 

Lop  Nor,  Sinkiang 

112 

252 

(Total) 

The  28,775  Station  magnitudes  have  been  fed  into  a  maximunvlikelihood  inversion  scheme  which 
simuttaneously  determines  the  event  size  and  the  station  correction,  as  well  as  the  specific  path  correc¬ 
tion  for  each  source-station  pair.  The  simultaneously-inferred  path  and  station  corrections  are  related  to 
known  geological/geophysical  features.  Applying  these  path  and  station  corrections  to  the  raw  station 
magnitudes  of  any  individual  explosion  yields  a  systematic  reduction  in  the  fluctuational  variation  of  sta¬ 
tion  magnitudes  across  the  whole  network  with  a  reduction  factor  ranging  from  1.2  to  3  for  all  Soviet 
events  in  our  data  set.  Most  Novaya  Zemlya  events  exhibit  a  variation  reduction  factor  of  2.  With  these 
path-corrected/station-corrected  mt,{Pmn)>  the  [NORSAR]  bias  between  the 

southwest  and  northeast  subregions  of  the  Soviet’;'  Batapan  test  site  is  assessed  as  0.07  magnitude  unit 
[m.u.],  which  is  significantly  smaller  than  that  of  previous  studies.  This  bias  can  be  further  reduced 
somewhat  when  the  mt,  based  on  the  first  motion,  mi,[Pa),  is  used.  First  motion  of  the  initial  short- 
period  P  waves  also  appears  to  be  a  very  favorable  source  measure  for  explosions  fired  in  hard  rock 
sites  underlain  by  a  stable  mantle  such  as  Semipalatinsk.  For  example,  based  on  nit,{Pa)  alone  and 
without  any  extra  cratering-to-contained  correction,  the  Balapan  explosbn  of  Jan  15, 1965,  is  estimated 
to  have  a  yield  of  119  kt.  The  mt,(Pa) -based  yield  estimate  for  the  JVE  event  of  Sep  14,  1988,  is  112 
kt.  Between  100  and  150  kt,  the  rCj^  bias  between  Eastern  Kazakh  and  NTS  using  our  mt,(Pmax)  values 
is  0.35  m.u. 
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2.  STATION  MAGNITUDE  COMPUTATION:  getmag 

“gaf/naj;"  computes  several  types  of  magnitudes  with  a  typical  command  of  the  following  form: 

gaunag  -IPhasa]  {-a  AmplUudeJ  l-p  Periodt  [-gh  Phase]  l-o  Origk^  /-•  Staton] 

All  the  arguments  required  are  read  through  the  command  line.  The  arguments  include  the  displacement 
amplitude  (-a)  in  nm,  the  period  (-p)  in  seconds,  the  phase  name  (-ph)  (e.0. ,  nit,.  Mg.  L^,  PS),  the  ori¬ 
gin  information  (-o)  which  includes  the  epicenter  and  the  event  name.  Each  phase  has  a  specific  formula 
for  determining  the  magnitude,  and  hence  different  argumerts  might  be  required.  The  formulae  are 
described  briefly  in  the  following: 

[1]  mt,  ■  log(A/T)  -t-  B(A)  for  20**  <  A  <  95°,  where  B(A)  is  the  (^stance  normalizer  derived  by  Veith  and 
Clawson  (1972). 

[2]  nit,{P„)m  log(A)  +  2.42  log(A)  -  3.95  for  A  <  10°  (cf.  Vergino  and  Mensing,  1990). 

[3]  P-wave  spectral  magnitude,  PS  ■  log(A)  +  0.5  iog[tan(lo)^(A)]  -t-  0.5  log  [d(lo)/d(A)]  for  20°  <  A  < 
100°  {cf.  Bullen  and  Bolt,  1985).  The  take-off  angle,  lo.  is  approximated  by  a  fourth  order  polynomial  in 
A  {cf.  Rivers  et  ai,  1980) 

[4]  For  Mg ,  two  different  formulae  are  used: 

If  A  >  25°  .  Mg  =  log(A^  +  1.66  log(A)  +  3.30  {cf.  lASPEI,  1967). 

If  10°  <  A  <  25°,  Mg  m  log(A/T)  +  1.07  log(A)  +  4.16  {cf.  NuttB  and  Kim,  1975). 

[5]  For  mt,(Lp) ,  Jih  and  Lynnes  (1993)  suggest  the  following  fonnula: 

.  4.0272  4  lo9A(4)  4  ilogW  4  |lofl(sln(  ^  )l  4  .  |1) 

Although  it  might  appear  to  be  different  from  most  other  formulae  in  use,  this  equation  is  actually 
equivalent  to  Nuttli's  (1986ab,  1987)  and  it  is  more  convenient  to  use.  For  instance,  a  seismic  source 
with  1-sec  Lg  amplitude  of  110  pm  at  10  km  epicentral  distance  would  correspond  to  a  nii,{Lg)  of 
4.0272  +  2.0414  +  0.3333  -  1.4019  +  0.0000  »  5.000,  the  same  value  that  Nuttli’s  original  2-step  formu¬ 
lae  would  give.  The  Qq  and  q  values  built  into  the  code  “getmag"  are  listed  in  Section  6. 

Example 

Sample  calls  of  “getmag"  such  as 

getmag  -mb  -a  7.3  -o  60.0  78.8  Event_1  -s  QUA  -p  0.9  -ph  Pa 
getmag  -Ms  -a  400  -s  BKS  -p  20.0  -o  37.0  - 170.0  Event_2  -ph  LR 
getmag  -PS  -a  100  -o  37 -10  Event_3  -s  TUC  -p  21.3  -ph  PSPE  -x  O.S-2.0 
getmag  -Lg  -a  0.3  -v  3.S  -o  50.0  78.8  Event_4  -s  KON  -p  0.9  -ph  CLg  -n  0.0 

Should  give 

mb(Pa)=  4.674  -o  60.000  78.800  Event J  -a  7.3  -p  0.90  -s  QUA  -ph  Pa 

Ms(LR)=  4.216  -c  37.000  -170.000  EventJ  -a  400.0  -p  20.0  -s  BKS  -ph  LR 

PS(PSPE)=  6.344  -o  37.000  -10.000  EventJ  -a  100.000000  -s  TUC  -ph  PSPE  -x  O.S-2.0 

mb(CLg)=  4.201  -o  50.000  7B.mO  EventJ  -a  0.3  -p  0.90  -v  3.50  -s  KON  -ph  CLg  -QO  700  -eta  0.40 
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3.  SINGLE-EVENT  MAXIMUM-LIKELIHOOD  ESTIMATOR:  emils  (domie) 

The  problem  of  estimating  body-wave  magnitudes  {mt, )  using  amplitudes  read  at  a  number  of 
recording  stations  is  frequently  complicated  by  the  fact  that  the  data  may  be  heavily  censored.  This 
arises  either  because  of  clipping,  where  all  amplitudes  can  be  determined  only  to  exceed  a  given  lower 
bound  (i.e.  the  right-censored  case  in  statistical  terms),  or  because  the  signals  are  weaker  than  the 
ambient  noise  level  and  hence  are  not  detected  (/.e.  the  left-censored  case).  If  one  simply  averages  the 
magnitudes  for  those  stations  which  detected  an  event,  without  regard  for  those  that  clipped  or  did  not 
record,  serious  biases  may  result  in  the  event  magnitude  estimated. 

For  single-event  network  mt,  determination,  at  least  three  types  of  station  magnitude  ought  to  be 
considered: 

(0]  the  station  magnitude,  X,  is  known  as  Xo, 

[1]  X  is  only  known  to  be  less  than  certain  level,  say,  tt, 

[2]  X  is  only  known  to  be  larger  than  certain  level,  say,  tz- 

We  assume  that  the  observed  station  magnitude,  X,  can  be  represented  as  the  sum  of  the  unk¬ 
nown  event  magnitude,  p,  and  a  perturbing  random  noise,  v, 

X  =  4  +  V  (21 

where  v  is  assumed  to  be  a  Gaussian  random  variable  with  mean  zero  and  standard  deviation  o. 
Elegant  maximum-likelihood  theory  can  be  derived  for  this  linear  model.  Suppose  there  are  no,  n^,  and 
n2  station  recordings  for  each  type,  respectively-  The  conditional  likelihood  function  of  the  censored 
observations  (  Xq,  tt,  ta)  given  the  network  magnitude  p  and  a  is 

L  (  Xo,  ti,  ta  I  p.  a )  =nP(  X|  =  Xoj  I  p,  a  )  *  nP(  <  tij  I  p,  o  )  *  nP(  Xj  >  taj  I  p,  a  )  ,  [3] 

i=i  j=i  i=i 

and  the  log-likelihood  function  is 

In  L  (  Xo,  tt.  tz  I  p,  a  )  =  -  -^  ln(2rto2)  -  ^ (^oj "  4)^  +  Z In  «I>(Ztj)  +  Zin  <I>(-Zai)  ,  [4] 

*  2<r  j  a  1  j  =  t  j  =  t 

where  Z|  s  (t|  -  p  )/o  ;  Xq,  tt,  and  ta  are  collections  of  the  observed  station  magnitudes  of  each  type, 
respectively,  and 

•Mu)  s  ^exp(^) ,  <l>(u)  =  £(Mx)dx  .  [5] 

are  the  probability  density  function  and  probability  distribution  function,  respectively,  of  the  standard  nor¬ 
mal  random  variable. 

Solving  s  0  implies  that  d,  the  optimal  estimate  of  a,  must  satisfy  the  following  necessary 
da 

condition; 
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6^  = 


Z{xoi  -  A  )* 


•Za 


dInL 


Solving  a  o  implies  that  A.  the  optimal  estimate  of  p.  must  satisfy  the  following  necessary 
dp 


condition: 


n-a-?,  *av 


Adding  (n^  +  n2)  A  to  both  sides  of  [6],  and  then  di^^ing  both  sides  by  (no  +  n^  +  n2)  yields 

1 


A  = 


rio  +  Hi  +  n2 


H  pi  ®\Zij;  j.1  ®(-Z2j) 


m 


18] 


The  right-hand  side  of  Equation  [7]  happens  to  be  the  sample  mean  of  “all”  data  with  the  censored 
measurements  replaced  by  their  corresponding  best  fill-in  (see  Appendix): 


—— (XE(XlX  =  Xoi]  +  £E[XlX<t,jl-H2EIXlX>t2i]).  [9] 

no  +  n,  +  nj  j.,  j., 

Consequently,  within  the  context  of  Gaussian  assurr^ion,  one  can  translate  those  seemingly  not-that- 
precise  statements  of  X  >  t  or  X  <  t  into  quantitative  constraints  which  can  couple  with  other  measure¬ 
ments  of  type  0  easily.  Thus  Equations  [8]  and  [9]  provide  the  theoretical  justification  of  an  iteration  pro¬ 
cedure  to  be  discussed  below. 


An  iterative  procedure  called  “EM  algorithm”  [Expectation-Maximization  algorithm]  (Dempster 
et  aL ,  1977)  can  be  applied  to  solve  for  p  and  a  in  a  very  straightforward  manner.  To  start  the  iteration, 
one  needs  an  initial  guess  of  a  and  p.  A  good  initial  value  of  p  is  the  sanple  mean  of  all  type-0  station 
magnitudes.  Since  bulletin  m^,  typically  exhibits  a  a  (of  single  obsenration)  around  0.3  magnitude  unit, 
this  value  can  serve  as  the  initial  value  of  a.  The  iteration  procedure  follows: 

[1]  Based  on  the  current  estimates  of  p  and  o,  replace  all  the  censored  data  with  their  corresponding 
conditional  expectations  {cf.  the  right-hand  side  of  Equation  [7]).  This  is  the  so-called  "E  step”  of 
the  EM  algorithm. 

[2]  Compute  A  as  the  sample  mean  of  these  “refined  observations". 

[3]  Update  the  estimate  of  a  using  Equations  [6]. 

[4]  Repeat  [1]-[3]  until  some  convergence  criterion  is  met. 

Steps  [2]  and  [3]  constitute  the  “M  step”  of  the  EM  algorithm.  Note  that  In  the  non-censoring 
case,  i.e. ,  n,  =  nz  =  0,  A  and  d  would  reduce  to  the  regular  sanple  mean  and  the  RMS  residual,  respec¬ 
tively: 
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[10] 


no  no 

2>o,  Bxoj-H)^ 

H - 

no  no 


Example 


The  algorithm  descr&ied  above  has  been  imptemented  as  a  utility  program  “en^s”  (“domie")  whidi 
expects  to  read  just  two  columns  of  data  representing  the  data  type  or  “>”)  and  the  actual 

data.  Take  Novaya  Zemlya  event  66300  (October  27,  1966)  as  an  example.  Table  2  lists  the  station 
mt,{Lg)  values  of  this  Novaya  Zemlya  event  based  on  our  mo(Lp)  formula  (cf .  Section  2)  as  well  as  the 
path  corrections  we  installed  {cf.  Section  6). 


Table  2.  Station  Recordings  of  Novaya  Zemlya  Explosion  66300 


Station 

A" 

Amplitude  [nm] 

COP 

24.57 

870.5 

KEV 

9.48 

<1833.6 

NUR 

17.22 

867.5 

STU 

31.67 

234.3 

UME 

15.58 

1168.2 

ESK 

29.23 

155.7 

1ST 

34.70 

49.5 

KON 

21.91 

789.3 

TRI 

33.38 

163.1 

Period  [sec] 


1.21 


0.88 


Velocity 


3.5 


mb(Lg) 


6.341 


<6.506 


6.389 


6.514 


6.525 


6.423 


6.464 


6.518 


417  0.24  3.6  6.221 


Amplituda  measurenwnt*  (umifhed  by  Fttvart  at  aL  (1993). 

There  are  8  good  signals  and  1  noisy  measurement: 

6.341 
"<"  6.S06 
”="  6.389 
6.514 
6.525 
6.423 
6.464 
6.518 
"="  6.221 

If  the  censored  recording  of  6.506  at  the  station  KEV  is  ignored,  the  event  magnitude  would  be 
6.42410.037.  The  program  “emils"  gives  the  maximum-likelihood  estimate  as  6.42010.034,  using  all  9 
observations.  Basically,  what  the  maximum-likelihood  method  does  is  to  utilize  the  censored  information 
of  mt,{Lg){KEV)  <  6.506  as  an  extra  constraint  to  refine  the  inferred  parameter  obtained  with  the  stan¬ 
dard  feast  squares.  For  this  event,  Nuttli  (1988)  gave  a  mi,(Lg)  of  6.45. 
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4.  SIMULTANEOUS  INVERSION  OF  EVENT  m^.  STATION  TERMS,  AND  PATH  CORRECTIONS 

4.1  General  Concepts  of  Joint  Inversion  Model 

As  described  in  Section  2,  the  conventional  definition  of  the  station  magnitude  is  computed  as 

irib  =  logio{A/T)  +  B(A)  .  [11] 

where  A  is  the  displacement  amplitude  (in  nm)  and  T  is  the  predominant  period  (in  sec)  of  the  P  wave. 
The  B(A)  is  the  distance-correction  term  that  compensates  for  the  change  of  P-wave  amplitudes  with 
distance  (e.g. ,  Gutenberg  and  Richter,  1956:  Veith  and  Clawson,  1972).  in  [11]  is  also  denoted  as 
mi  in  Marshall  et  al.  (1979).  The  ISC  bulletin  /r%  is  just  the  network  average  of  these  raw  station  m^ 
values  without  any  further  adjustment.  That  is,  we  assume  a  linear  model  as  the  following: 

^i>(D  =  E  +  v(j)  ,  [12] 

where  m^  (j)  is  the  station  magnitude  recorded  at  the  station  j  for  the  event  of  size  E,  and  v  is  the  ran¬ 
dom  perturbing  term. 

Now  consider  Ne  explosions  detonated  at  Nf  source  regions  that  are  recorded  at  some  or  all  of  Ns 
stations.  In  LSMF  [Least  squares  Matrix  Factorization]  and  the  standard  GLM  [General  Linear  Model] 
schemes  {e.g.,  Douglas,  1966;  Blandford  and  Shumway,  1982;  MarshaU  et  al.,  1984;  LilwaB  et  ai, 
1988;  Jih  and  Shumway,  1989;  Murphy  et  aJ.,  1989),  it  is  assumed  that  the  observed  station  m^,  (i,j)  is 
the  sum  of  the  true  source  size  of  the  Mh  event,  E(i).  the  receiver  term  of  the  j-th  station,  SO),  and  the 
random  noise,  v(i,j); 

ff»(,(i.j)  =  E(i)  +  S(0  +  v(i,j)  ,  [13] 

The  receiver  term,  S(j),  is  constant  with  respect  to  all  explosions  from  different  test  sites,  and  hence  it 
would  inherently  reflect  the  “averaged"  receiver  effect  ~  provided  the  paths  reaching  the  station  have 
broad  azimuthal  coverage.  When  world-wide  explosions  are  used,  the  standard  deviation  (o)  of  the 
noise  v  in  [13]  is  typically  around  0.3  m.u. 

if  LSMF  or  GLM  is  applied  to  events  within  a  smaller  area  of  source  region,  then  the  a  of  v  in  [13] 
could  reduce  to  0.15-0.2  m.u.  However,  the  result  of  such  “single-test-site  GLM"  approach  should  be 
interpreted  or  utilized  cautiously.  The  event  mt,  values  (/.e.,  the  “E"  term  in  [13])  so  determined  are 
excellent  estimates  of  the  “relative  source  size"  for  that  test  site  only.  If  this  “single-test-site  GLM" 
inversion  is  applied  to  several  test  sites  separately,  ft  may  not  be  easy  or  obvious  to  find  a  consistent 
baseline  for  estimating  the  “absolute  yield",  since  the  recording  network  is  typically  different  from  one 
test  site  to  another,  and  hence  the  station  terms  are  inevitably  inconsistent.  Furthermore,  the  station 
terms  derived  by  the  “single-test-site  GLM"  may  not  necessarily  represent  the  attenuation  underneath 
the  receiver  side  alone.  They  could  be  “contaminated"  or  sometimes  even  overwhelmed  by  the 
path/near-source  effects  shared  by  the  explosions  confined  in  a  narrow  azimuthal  range.  This  could 
explain  the  once  puzzling  and  controversial  phenomenon  Butler  and  Ruff  (1980)  (also  Butler,  1981;  Bur¬ 
dick,  1981)  reported,  namely  that  using  Soviet  explosions  from  one  test  site  alone  may  fail  to  discern  the 
attenuation  differential  between  the  eastern  and  western  U.S.  There  is  no  doubt,  however,  that  the  GLM 
or  LSMF  type  of  methodology  can  infer  the  station  terms  which  are  strongly  correlated  with  the  upper 
mantle  attenuation  underneath  the  stations,  provided  the  seismic  sources  have  a  broad  spatial  coverage 


7 


as  did  those  in  North  (1977),  Douglas  and  Marshall  (1983),  Lilwall  and  Neary  (1985),  Ringdal  (1986),  Jih 
and  Wagner  (1991),  and  many  others.  The  event  magnitude  derived  with  Equation  [3]  is  hereby  denoted 
as  rTf22-  In  Marshall  et  at.  (1979),  a  priori  information  about  the  P„  velocity  underneath  each  station  is 
used  to  determine  its  associated  “deterministic"  receiver  correction,  S(j),  and  the  network-averaged 
magnitude  based  on  the  station-corrected  magnitudes  is  called  m2  ■  The  receiver  corrections  as  derived 
in  Equation  [3],  however,  are  inferred  jointly  from  a  suite  of  event-station  pairs,  and  no  a  priori  geophy¬ 
sical  or  geological  condition  is  assumed  (and  hence  the  different  notation  m22)-  The  high  correlation 
between  the  tectonic  type  and  the  GLM  station  terms  suggests  that  the  empirical  station  corrections  do 
reflect  the  averaged  upper  mantle  conditions  underneath  the  receivers,  if  the  azimuthal  coverage  at  each 
station  is  broad  enough. 

Jih  and  Wagner  (I992ab)  propose  to  reformulate  the  whole  rrxxJel  [13]  as 

/n*  (i.j)  =  E(i)  +  SG)  +  F(k,j)  +  v(i,j)  ,  [14] 

where  F(k,j)  is  the  correction  term  at  the  j-th  station  for  the  propagation  effect  or  the  near-source 
focusing/defocusing  effect,  which  is  constant  for  all  events  (including  this  ith  event)  in  the  k-th  “geobgi- 
cally  and  geophysically  uniform  region".  For  each  seismic  statbn,  this  F  can  be  regarded  as  its  azimu¬ 
thal  variation  around  the  mean  station  term  S.  However,  as  explained  prevbusly,  it  would  be  more 
appropriate  to  consider  F  the  path  or  near-source  term  because  the  back  azimuths  at  the  station  could 
be  nearly  identical  for  adjacent  test  sites  (such  as  Degelen  and  Murzhik),  and  yet  the  “F"  terms  could  be 
very  different.  By  incorporating  t^"  F  term  into  the  model,  the  o  for  world-wide  explosions  is  reduced  to 
about  0.2,  roughly  the  same  level  that  which  a  “single-test-site  GLM”  could  achieve.  Intuitively,  the 
present  scheme  (Equation  [14])  provides  a  more  detailed  (and  hence  better)  model  than  that  of  Equatbn 
[4]  in  describing  the  whole  propagation  path  from  the  source  towards  the  receiver.  Simply  put,  Equatbn 
[13]  yields  a  stronger  fluctuatbn  in  the  source  terms,  E,  as  well  as  a  larger  standard  deviatbn  of  v 
because  each  term  in  the  right-hand  side  of  Equation  [13]  woub  have  to  "absorb”  part  of  the  missing  F 
term.  The  resulting  new  event  magnitude  (viz.,  E(i)  in  [14])  is  hereby  called  frizi  to  avoid  confusbn 
with  the  rriz  defined  in  Marshall  et  al.  (1979)  that  corrects  for  the  source-regbn  attenuation  and  statbn 
terms  solely  based  on  published  P„  velocity. 

Roughly  speaking,  the  model  described  in  [14]  has  the  foibwing  advantages; 

•  It  provbes  more  stable  mi,  measurements  across  the  whole  recording  network,  as  compared  to 
the  conventional  GLM  or  LSMF  procedure  whbh  only  corrects  for  the  statbn  terms.  The  reductbn 
in  the  standard  deviatbn  of  network  mi,  from  rrix  to  Wzs  could  reach  a  factor  of  nearly  3.  As  a 
result,  the  scatter  in  nizs  versus  log(yield)  is  smaller  than  that  for  other  m^ . 

•  The  separatbn  of  the  path  effect  from  the  statbn  effect  is  a  crucial  step  to  investigate  the  varbus 
propagation  phenomena,  whbh  in  turn  woub  improve  our  urberstanding  of  the  seismb  source  as 
well. 

We  have  applied  this  model  to  252  worldwide  expbsions,  and  the  resulting  rriza  values  of  these 
explosions  are  listed  in  Table  3.  The  132  statbns  are  selected  such  that  each  station  records  10  or 
more  good  explosion  signals.  There  are  relatively  fewer  explosions  recorded  at  the  modem  digital 
stations/networks.  As  a  result,  WWSSN  is  still  the  core  recording  network.  In  this  data  set,  there  are 
16,716  signals,  10,055  noise  measurements,  and  2,004  clipped  measurements  from  18  test  sites  that 


8 


are  used  to  invert  for  the  3,269  unknown  parameters  with  the  maximum-likelihood  approach.  The  stan¬ 
dard  deviation  of  v(i,j)  in  [14]  is  0.189,  as  compared  to  that  of  0.281  if  the  conventional  GLM  (Equation 
[13])  is  applied  to  the  same  data  set.  The  algorithm  and  sample  input  files  are  described  in  the  next  sec¬ 
tion. 

4.2  Maximum-llkellhood  General  Linear  Model:  migim 

"migim"  simultaneously  inverts  for  the  maximum-likelihood  estimate  of  event  magnitudes  and  sta¬ 
tion  corrections,  as  well  as  the  path  terms  with  a  data  set  of  which  some  stations  might  fail  to  detect  the 
signal  (due  to  the  noise  contamination)  or  might  be  clipped  due  to  the  limited  dynamic  range.  It 
assumes  a  general  linear  model  [GLM]  of  the  form: 

X(i,j)  =  E(0  +  SO-)  +  F(k,j)-fv(i.|)  ,  [15] 

where  E(i),  S(j),  and  F(k,j)  are  the  unknown  source  size  of  the  i-th  event,  the  station  term  at  the  j-th  sta¬ 
tion,  and  the  propagation  effect  from  k-th  test  site  (at  which  the  i-th  event  is  located)  to  the  j-th  station, 
respectively-  X  (i,j)  is  the  observed  station  magnitude  of  event  i  as  observed  at  the  station  j.  The  pro¬ 
gram  also  has  an  option  to  solve  for  E(i)  and  SO)  for  a  simpler  linear  model: 

X(ij)  =  E(i)-KS(j)-hv(i,j)  .  [16] 

Here  we  assume  that  there  are  four  types  of  data  available: 

[0]  the  observed  magnitude,  X,  is  known  as  Xo, 

[1]  X  is  only  known  to  be  less  than  certain  level, 

[2]  X  is  only  known  to  be  larger  than  certain  level,  and 

[3]  X  is  missing. 

As  discussed  in  the  previous  section,  the  resulting  event  magnitude,  E(i),  in  Equations  [15]  and 
[16]  is  called  /n2.g  and  m2.2.  respectively. 

Sample  Input  Files 

"migim"  reads  in  a  file  “Events"  which  specifies  the  input  parameters  as  well  as  all  the  event  files 
that  will  be  required  in  the  GLM  inversion. 

#  1—  give  output  file  name 

GLMjest 

»2—  give  label 

Testing  sager  geole^bin/s^mlglm,  OS  13 93 

»3~-  estimator  [O.SsLSMF,  1,4=ULE  (recommended),  2,S=ILS] 

1 

S4~-  how  many  events  should  a  good  station  record?  (1,  2,  3, ...) 

2 

*5—  give  distance  Hag  S  acceptable  distarxe  window 
1  20.00  95.000 

#6—  choose  terse  level  (0, 1,2,...) 

1 

83230.1 60958.nz.pmax  NZB30B18  73.38n  54.91e  NNZ 
83268. 130957.nz.pmax  NZ83092S  73.35n  54.50e  NNZ 
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a42a0.0629S7jtz.f>max  NZB4102S  73.3/n  54.969  NNZ 
e7214.01S9S»JU.pmm»  H26706O2  73.S4n  54.629  NNZ 
88l2e.224958.m4m9X  NZ860S07  73.36n  54.449  NNZ 
88338.051 9S3.nz.pmax  NZ881204  73.38n  SS.OOa  NNZ 

The  first  portion  of  the  input  parameter  fiie  is  seif-expianatoiy.  There  are  aclualiy  6  e^imators  to 
choose  from: 

0  MoMng  ITI22  *»»8»  LSUF  (Bquaaho  (16^ 

1  soMng  ni22  ^  fEquaion  [16]) 

2  soMng  ¥irilh  ILS  (EqiMion  [1^ 

3  soMng  fn2»  M61  LSiJF  (Equnion  [15J) 

4  soMng  fflza  Mth  ULE  (BqunHon  [15]) 

5  soMng  1712.2  Mth  kS  (Equation  [15J) 

Estimators  0  through  2  are  suitabie  for  the  case  of  a  single  test  site,  or  if  the  path  effects  from 
different  test  sites  are  to  be  ignored.  Estimators  3  through  5  are  suitsdsle  for  the  case  of  muttiple  test 
sites.  The  maximum-likelihood  estimator  [MLE]  has  received  more  attention  than  has  the  aRemate 
method  of  the  iterative  least  squares  [ILS].  The  methodobgicai  similarities  and  differences  between 
these  two  methods  are  discussed  in  detaH  in  Jih  arxf  Shumway  (1989). 

Each  of  the  remaining  fines  in  the  input  file  specifies  the  event  file  name,  date,  event  name,  geo¬ 
detic  coordinate  of  the  event,  and  the  test  site  with  free  format  (no  quotation  marks!). 

Each  event  file  {e.g. ,  83230.160958.nz.pmax)  contains  a  list  of  stations  as  wefi  as  the  correspond¬ 
ing  measurement  with  format  (a6,a1  ,f5.3)  as  shown  in  the  following  sanple  file: 

ALO  >5.113 
ANMO  <5.274 
BJI  5.411 
BKS  S.6B8 
BOD  6.040 

KM!  9.555  «—  to  bo  rejected 

Any  fields  after  the  12th  byte  are  generally  ignored  except  when  the  14th  byte  is  a  '#'  sign.  In  that 
case,  this  record  will  be  totally  rejected.  This  feature  is  especially  useful  when  quality  control  is  knposed 
on  of  the  input  data. 

The  routine  also  needs  a  listing  of  stations  (called  ’List’)  in  the  free  format  (thus  the  station  codes 
must  be  in  quotes!).  Only  3  columns  are  needed.  GLM  will  stop  and  remind  the  user  if  the  coordinate 
of  a  station  is  missing  or  if  some  event  has  no  signal  at  all. 


"AAE- 

9.0291660 

38.765556 

"AAM- 

42.299721 

-83.656113 

“AKU" 

65.686666 

-18.106667 

"ALQ" 

34.942501 

-106.457497 

"ANMO" 

34.946194 

-106.456665 

"ANTO" 

39.900002 

32.783333 
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4.3  Iteration  Procedure 


Equations  [15]  and  [16]  are  special  cases  of  general  linear  models  [GLM].  An  Kerative  procedure 
based  on  the  EM  algorithm  is  presented  below.  The  basic  ideas  are  very  sintilar  to  those  underlying  the 
single-event  network  averaging  presented  in  Section  3. 

Step  0 

Set  up  initial  conditions  as  follows: 

[1]  a  z  0.3  magnitude  unit, 

[2]  SO)  H  0forj  =  1.2 . Ns, 

[3]  F(k,j)  3  0  for  i  =  1,  2 . Ns,  and  k  =  1,  2 . Np, 

Step  1 

Compute  event  magnitudes,  E(i),  for  i  -  1 . Np  as 

E(0  =  -i|jj-BX(i,j)  -  SO)  -  F(k,j)] , 

where  #0)  is  the  number  of  stations  that  “recorded”  the  event  i. 

Step  2 

Compute  station  corrections,  SO),  for  j  -  1 . Ns  as 

S0  =  ’ 

where  «(i)  is  the  number  of  events  “recorded”  at  station  j. 
steps 

Compute  path  corrections,  F(k,|),  for  j «  1 . Ns;  k  >  1,...,  Np  as 

F(k.D  =  -  E{i)  -  SO)] . 


where  #(  (k,j) )  is  the  number  of  paths  from  the  test  site  k  (where  the  event  i  is  located)  to  the  sta¬ 
tion  j.  This  step  is  skipped  if  options  0  through  2  are  chosen.  Consequently,  F(k,  j)  will  remain  0 
for  all  k  and  j  when  /n2^  is  the  desired  event  magnitude. 

Step  4 

Remove  the  mean  of  SO)  from  each  station  term  so  that  £  SO)  =  0  . 

i 

Steps 

For  each  source-station  pair,  (i,  j),  compute  p(i,j)  s  E(i)  +  SO)  +  F{k,j). 

Step  6 

For  estimators  1  and  4,  compute  o(MLE)  via 


"o 

-  Moj 

_ H _ 

_  <Mzij)  ^  'Wzzi) 


For  estimators  2  and  5,  compute  o(ILS)  via 
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[181 


a2  = 


£(xoj  -  m<k)^ 
H _ 


^Zlj) 


♦(ZZi) 


no  +  n, +  02  ^''*K-Z2i)' 

where  Z)  s  (ti  -  ^  )/a  ;  Xq.  ti,  and  t2  are  collections  of  the  observed  station  maonitudes  of  each 
type,  respectively,  and 


are  the  probability  density  function  and  probability  (fistribution  function,  respectively,  of  the  stan¬ 
dard  normal  random  variable. 

For  estimators  0  and  3,  [17]  and  [18]  would  be  equally  applicable  since  n^  -  n2  »  0. 

Step? 

Replace  censored  and  missing  obsenrations  X(i,j)  with  the  corresponding  conditional  expectations: 
For  type-1  paths:  E  ( X  I  X  <  t,j  1  =  p  -  a  . 

WtZij) 


For  type-2  paths:  E[XlX>tql  =  p  +  a  ^  . 

®V~Z2jJ 

For  type-3  paths:  E  [  X  I  X  is  missirtg  1  ■  p(l,j)  =  E(l)  +  S(j)  +  F(k,j)  . 

These  conditional  expectations  are  then  used  as  X(i,i)  in  steps  1  through  3. 
steps 

Repeat  steps  [1]-[7]  to  update  E,  S,  F,  and  o  until  convergence. 

In  the  first  iteration,  only  type  0  data  are  used  In  steps  1  through  3.  Starting  from  the  second  loop, 
however,  all  types  of  observations  are  used  with  censored  data  replaced  by  their  corresponding  “refined 
pseudo-observations"  as  des<;.ibed  in  step  6.  In  other  words,  the  symbol  X(i,j)  in  steps  1-3  actually 
represents  the  condKional  expectation  of  X  given  the  censoring  or  non-censoring  assumption.  For  type- 
0  data,  E  [  X  I X  -  xoj  ]  -  Xo|,  and  hence  the  actually  obsenred  magnitude  is  utilized  in  each  iteration 
without  change.  For  other  types  of  data,  however,  ttw  “expected”  observation  will  be  varying  as  the 
iterations  proceeds,  since  the  optimal  estimate  of  cr  and  all  other  parameters  will  change  at  each  step. 

Once  the  “E  step"  {viz. ,  steps  5  and  7)  is  executed,  the  “M  step"  {viz. ,  steps  1  through  4)  in  each 
iteration  loop  can  be  replaced  with  standard  matrix  inversion  techniques  such  as  Singula'  Value  Decom¬ 
position.  [SVO]  or  Gaussian  elimination  method.  To  do  so,  type-3  paths  should  be  excluded  from  step  7. 
Numerical  algorithms  like  SVO  and  Gaussian  elimination  are  called  direct  methods.  However,  direct 
methods  can  be  impractical  if  the  design  matrix  is  large  and  sparse.  In  our  case,  the  linear  system 
involves  3,269  unknown  parameters  and  28,775  station  magnitudes.  For  these  types  of  problems,  itera¬ 
tive  methods  are  superior  to  Gaussian  elimination  and  matrix  factorization.  The  largest  area  for  the 
application  of  iterative  methods  is  that  of  the  linear  systems  arising  in  the  numerical  solution  of  partial 
differential  equations.  Systems  of  orders  10,000  to  100,000  are  not  unusual  in  aerospace  sciences, 
although  the  majority  of  the  coefficients  of  the  systems  are  typically  zeros. 
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4.4  Examplo:  Running  GLM  with  Geotech’s  Whole  Dataset 

The  following  script  runs  GLM  with  our  whole  data  sM.  Since  we  adopted  a  very  restrictive 
criterion  in  screening  the  stations,  some  events  could  have  no  good  signal  after  the  quality  check.  The 
program  “migim"  returns  a  message  indicating  ttiat,  and  it  slope.  The  user  can  then  manualy  edit  the 
listing  of  events  to  delete  such  events  and  then  resubmit  the  GLM  job.  The  following  script  mckJdes  a 
section  to  perform  this  function  sujtomaticaay,  based  on  R.  R.  Baumstark's  suggestion. 

*—  script  to  run  GLM  kworsion  mth  TG's  whole  date  set:  nas  needed:  “UsT.  “£VJIsr. ... 

LOOP: 

cat«  EOP>  Eventsjtead 

*1.  Give  output  He  name  (a13)  (wil  overwrite  "GLMjnsgT) 

GLMjout 

*2.  Give  a  label  (aSO) 

WWSSN  mb  inversion  with  GLM 

53.  Choose  estimator  (LSMF:0,3:  MLE:1,4:  ILS:2,S) 

4 

54.  At  least  how  many  avertts  should  a  good  station  record?  (t,2,3....) 

10 

#5.  Distartoe  Hag  (1:  on,  0:  ofO  4  min.  max.  distatKe  (in  deg.)  acceptable 
1  20.0  95.000 

S6.  Choose  terse  level  of  output  (0, 1,2,...) 

1 

EOF 

cat  Eventsjiead  EV_Sst>  Events 

S—  run  mIgIm:  if  execution  not  complete,  error  message  would  stit  be  'GLMjns^, 

5  inrScating  some  user  intervention  might  be  needed. 
mlglm3 

•—  QC  loop  (added  May  10,  1993,  based  on  Boottter’s  suggmion) 
s  to  delete  "bad"  events  tiom  EVJisL  and  re-run  "mlglnf 
if  ( -a  GLMjnsg )  then 

iff  'grep  ’has  no  signal'  GlM_msg\wc  -I  *>  0)then 
grep  has  no  sigrtaT  GLMjnsg  >  FOO 

echo  Rejecting  "wc  -I  FOOlawk  '{print  $1}"  events  with  no  signals: 
cp  EVJist  Keep 

foreach  bad  ( 'avrit  Iprint  $6] '  FOO" ) 
echo  '  'rejecting  event  $bad 
awk  f  if($1  l‘‘'"Sbad''  )  print  $0]'  Keep  >  foo 
mv  foo  Keep 
end 

mv  Keep  EVJist ;  rm  FOO  GLMjnsg 
ifCwc-l  EV_list\awk  {print  $1]"  t^O)  then 
echo  NO  events  lefL  exiting, 
exit 

endif 

echo  fterunning  GLM  on  reduced  data  set 
goto  LOOP 

endif 

endif 

S—  end  of  QC  hop 
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4.5  Tables  of  Resulting  Event  m,,  Values 

The  “migim'’  program  generates  six  ASCII  output  buffers  in  addition  to  the  error  message  file.  The 
following  script  reads  one  of  the  GLM  output  buffers,  ‘iort.48'\  and  makes  a  table  of  values  sorted 
by  test  sites.  The  buffers  “tome",  “tort39“,  and  “fort.49  the  resulting  station  and  path  corrections 
which  are  also  suitable  for  map  plottirtg  purposes  (cf.  Sections  4.6  and  4.7). 

»—  script  to  maka  a  tabla  at  mb  mluas 
set  nonomalch 

*—  separate  Pa,  Pb.  Pmax  bom  “nUgIm"  output  buffer 
sect  -n  -e  '1,24  Ip'  <  tort.48  >  GUd-Pa 
sad-n-e  '242.492p' <  hrt4a>  GtMPb 
sed-n-e  '403, 744p'  <  hrt4B>  QIMPmax 

#—  group  evartts  in  QIMPmax  by  site 

sed-n-e  1almen<bofjshoaftp'<  GLMPmas>  USA 

sed -n  -a  1azg22apt66IJpne20aug74lp' <  QIMPmax>  PNB 

sad -n  -a  1nnz2SoelS4IJsm10eet7Slp' <  QIMPmax*  NZ 

sad  -n  -a  1kon18dac66IJdak22mayBap'<  QIMPmax*  Deg*Mzk 

sad  -n  -a  IsaklSjanSSIJaekO^ulOOIp'  <  QIM.Pmax *  Batapan 

sed  -n  -a  1beryljch21may92lp'  <  GIM.Pmax  *  Other 

rm  GLM.Pmax 

if  ( -a  USTjnb )  rm  USTjnb 

set  siteH  USA  PNE  NZ  Other  Oegi^Mzk  Bafspmi ) 

foreachk(  1  23456) 

*—  for  each  test  site,  search  Pat  Pb  for  each  Pmax 
foreach  Hrte  (’  ‘  cat  $silel$IQ  "'  ) 
set  INm'  echo  $tne  ‘ 
set  IDm‘  echo  $IN[1] ' 

set  SNCEm'  echo  $lfH6]  $IN[7]  VNIB]  PN[3]  ’ 
set  mbs-'  echo  $IN[2]  ‘ 

( grep "StOr  GIMPa *  foundl  )*  6  fdevfnuU 
if  (-2  foundl  )  then 

set mb1^‘ echo' _ ‘ 

else 

sat  mb1=’  Boomer  foundl  2 ' 

erxSf 

( grep  'SIP'  GIM.Pb *  found2 }*  6  idedmiH 
if  ( -2  found2  )  then 

set  mb2=‘ echo ' _ ' 

else 

set  mb2='  Boomer  found2  2  ‘ 

endif 

echo  $10  SSNCE  Smbi  $mb2  $mb3  »  LIST  mb 
rm  found' 

end 

rm  $site[$k] 
end 

#---  end  of  loop  on  k 
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mOLIir 
caf«  />  DO 
«—  final  hmatting 

tnak  •iprini/  -  %$%3b  %3s  %3s%4Jlf%4^%4.^%4.2f  - _ )  *  <  USTjnb  >mbMI:m  LIST  mb 

I 

eah  DO ;  m  DO 

and 

and 

•—  and  of  hop  onmSn 
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Table  3.  Magnitudes  o<  Semipalatinsk  Explosions 


Event 


Date 


6501 15B 


6511210 


6602130 


6603200 


6605070 


6610190 


661218M 


6702260 


67091 6M 


670922M 


671122M 


680619B 


6809290 


690531 M 


6907230 


6909110 


691130B 


691228M 


700721 M 


701104M 


7103220 


7104250 


710606M 


710619M 


710630B 


711009M 


711 021 M 


7112300 


72021 OB 


7203280 


7208160 


720826M 


720902M 


721102B 


721 21 OB 


7212100 


730723B 


731214B 


#  of  Signals 


Ns  Nn 


461  2 


48  15  1 


51  4  10 


49  9  8 


9  26  1 


51  10  5 


55  8  1 


48  9  6 


36  29  2 


35  31  1 


7  64  0 


28  3  2 


50  8  6 


30  31  0 


38  21  1 


19  39  0 


50  0  0 


45  9  3 


Magnttudes  [rrizs] 


4314  3 


37  5  0 


38  12  2 


41  13  0 


31  19  1 


2712  3 


32  9 


1630 


34  8  2 


28  17  0 


23  23  1 


29  15  2 


1529 


44  2  11 


30  7  5 


531  1 


49  8  6 


0.03 

5.52 

5.75 

0.02 

4.98 

5.25 

0.02 

5.73 

5.98 

0.02 

5.43 

5.71 

0.03 

4.11 

4.25 

0.02 

5.43 

0.02 

5.42 

5.66 

0.02 

5.45 

5.70 

0.02 

4.69 

4.96 

0.02 

4.55 

4.87 

0.02 

4.01 

0.03 

4.72 

^^^2^211111111 

0.02 

5.24 

5.52 

0.02 

4.50 

4.91 

0.02 

4.73 

5.04 

0.03 

4.16 

4.40 

0.03 

5.41 

5.79 

0.03 

5.29 

5.58 

0.02 

4.72 

5.06 

0.02 

4.96 

5.17 

0.02 

5.13 

5.42 

0.03 

5.45 

5.71 

0.03 

4.91 

5.25 

0.03 

4.89 

5.19 

0.03 

4.37 

4.76 

0.03 

4.82 

5.05 

0.03 

4.91 

5.24 

0.04 

5.44 

0.03 

4.86 

5.12 

0.03 

4.50 

4.84 

0.03 

4.46 

4.75 

0.03 

4.72 

5.06 

0.03 

4.18 

4.44 

0.03 

5.62 

5.94 

0.03 

5.84 

0.03 

5.09 

5.41 

0.03 

5.76 

6.00 

0.02 

5.30 

5.59 

1 )  BSW  =  SW  subsite,  Balapan;  BNE  NE  subsite,  Balapan;  BTZ  =  transition  zone,  Balapan;  Deg 
zhik. 

2)  Ns  =  #  of  signals,  Nn  =  #  of  noise  measurements,  Nc  =  «  of  dips. 

3)  standard  error  in  the  mean. 


5.90 


5.46 


6.16 


5.93 


4.54 


5.61 


5.88 


5.93 


5.21 


5.15 


4.38 


5.30 


5.72 


5.14 


5.26 


4.72 


5.95 


5.78 


5.31 


5.38 


5. 


5.90 


5.45 


5.42 


5.04 


5.25 


5.47 


5.62 


5.35 


5.07 


5.00 


5.29 


4.71 


6.16 


6.03 


5.64 


6.18 


5.80 


Oegelen  Mountain;  Mzk  =  Mur- 


100-150 


29 


125 


20-150 


20-150 


20-150 


<20 


10 


<20 


<20 


<20 


16 


<20 


125 


46 


<20 


<20 


20-150 


90 


16 


<20 


<20 


12 


23 


20-150 


16 


6 


8 


<20 


2 


165 


140 


20-150 
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Table  3.  Magnitudes  of  Semipatatinsk  Explosions  (continued) 


Event 

#  of  Signals 

Magnitudes  (/n2.g] 

Date 

Site 

Ns  Nn  Nc 

S.E.M. 

P. 

Pb 

Pmn 

7504278 

8TZ 

18  1  1 

0.04 

4.90 

5.23 

5.53 

7607048 

8SW 

38  0  5 

0.03 

5.20 

5.82 

7611238 

8NE 

31  0  0 

0.03 

5.31 

5.68 

5.81 

7612078 

8SW 

172  1 

0.04 

4.96 

5.37 

5.60 

770329D 

Deg 

25  14  0 

0.03 

4.42 

4.80 

5.09 

7705298 

8SW 

30  4  4 

0.03 

5.25 

5.50 

5.71 

7706298 

8NE 

27  11  0 

0.03 

4.77 

4.94 

5.18 

770730D 

Deg 

21  16  0 

0.03 

4.31 

4.71 

4.96 

7709058 

8NE 

31  1  4 

0.03 

5.32 

5.57 

5.83 

780326D 

Deg 

25  6  0 

0.03 

5.01 

5.31 

5.54 

780422D 

Deg 

21  9  0 

0.04 

4.58 

4.84 

5.08 

7806118 

8SW 

17  0  1 

0.04 

5.26 

5.53 

5.83 

7807058 

8SW 

38  7  7 

0.03 

5.21 

5.48 

780728D 

Deg 

36  9  6 

0.03 

5.08 

5.38 

7808298 

8NE 

19  0  1 

004 

5.62 

5.90 

5.96 

7809158 

8SW 

37  1  6 

0.03 

5.44 

5.67 

5.84 

7811048 

8NE 

40  9  6 

0.03 

5.15 

5.39 

5.61 

7811298 

8SW 

30  0  1 

5.53 

5.83 

5.89 

7906238 

8SW 

40  3  3 

0.03 

5.65 

5.88 

6.08 

7907078 

8NE 

32  0  0 

0.03 

5.37 

5.63 

5.85 

7908048 

8SW 

40  5  20 

0.02 

5.89 

6.11 

8T2 

33  0  0 

0.03 

5.61 

5.90 

6.10 

7910288 

8NE 

44  5  13 

0.02 

5.51 

5.74 

5.97 

7912028 

8SW 

18  0  1 

0.04 

5.41 

5.67 

5.90 

7912238 

8SW 

41  3  17 

0.02 

5.60 

5.89 

6.13 

800522D 

Deg 

36  23  1 

0.02 

4.74 

4.99 

5.20 

8006298 

8SW 

46  6  6 

0.03 

5.21 

5.46 

5.67 

8009148 

8SW 

34  5  6 

0.03 

5.50 

5.83 

6.09 

8010128 

8NE 

27  0  0 

0.04 

5.75 

5.88 

8012148 

8TZ 

33  0  0 

0.03 

5.46 

5.75 

5.96 

8012278 

8NE 

29  0  0 

0.04 

5.56 

5.77 

5.92 

8104228 

8SW 

31  0  0 

0.03 

5.41 

5.70 

5.92 

8109138 

8TZ 

24  0  0 

0.04 

5.64 

5.93 

6.09 

8110188 

8SW 

41  4  7 

0.03 

5.50 

5.78 

5.99 

8111298 

8SW 

37  12  5 

0.03 

5.05 

5.32 

5.53 

8112278 

8SW 

29  0  1 

0.04 

5.67 

5.99 

6.18 

8204258 

8SW 

46  3  9 

0.03 

5.54 

5.80 

6.00 

8207048 

8SW 

25  1  0 

0.04 

5.68 

5.97 

6.08 
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Table  3.  Magnitudes  of  Semipalatinsk 


Event 


Explosions  (continued) 


Magnitudes  [/njg] 


820831 B 


821205B 


821226B 


83061 2B 


831006B 


831026B 


831120B 


840219B 


840329B 


840425B 


840526B 


840714B 


840915B 


841027B 


841202B 


841216B 


8412adB 


850210B 


85042SB 


850615B 


850630B 


850720B 


87031 2B 


870403B 


87041 7B 


870620B 


871115B 


871213B 


880403B 


880504B 


880614B 


8809 14B 


881112B 


881217B 


890708B 


Event 


(Date) 


NNZ250CT64 


NNZ270CT66 


NNZ210CT67 


NNZ07NOV68 


NNZ140CT69 


NNZ14OCT70 


NNZ27SEP71 


NNZ28AUG72 


NNZ12SEP73 


NNZ29AUG74 


NNZ23AUG75 


NNZ21CXJT75 


NNZ29SEP76 


NNZ20OCT76 


NNZ01SEP77 


NNZ09OCT77 


NNZ10AUG78 


NNZ27SEP78 


NNZ24SEP79 


NNZ18(XT79 


NNZ11OCT80 


NNZ01OCT81 


NNZ110CT82 


NNZ18AUG83 


NNZ25SEP83 


NNZ250CT84 


NNZ02AUG87 


NNZ07MAY88 


NNZ04DEC88 


NNZ24OCT90 


Table  3.  Magnitudes  of  Novaya  Zemlya 


#  of  Signals 


Explosions 


Magnitudes  [/njg] 


NsNnNc 


20  0  0 


56  0  13 


53  5  3 


58  1  5 


59  2  7 


35  0  22 


S.E.M. 


25  0  18 


27  0  12 


23  0  17 


27  4  7 


25  34  0 


25  2  2 


18  22  0 


39  3  18 


42  7  10 


33  2  16 


39  7  14 


42  4  6 


43  4  5 


32  11  5 


30  4  5 


31  4  5 


22  3  4 


24  3  6 


27  4  1 


20  4  2 


SNZ27SEP73 

48 

3 

1 

0.03 

SNZ270C73A 

14 

0 

24 

0.03 

SNZ270C73B 

9 

28 

0 

0.03 

SNZ270C73C 

4 

34 

0 

0.03 

SNZ02NOV74 

12 

0 

29 

0.03 

SNZ180CT75 

21 

0 

21 

0.03 

0.04 

4.43 

4.51 

0.02 

6.07 

6.30 

0.02 

5.42 

5.61 

0.02 

5.60 

5.84 

0.02 

5.76 

5.96 

0.03 

6.43 

6.62 

0.03 

6.24 

6.43 

0.03 

5.98 

6.23 

0.03 

6.36 

6.67 

0.03 

6.13 

6.38 

0.03 

6.15 

6.38 

0.03 

6.10 

6.33 

0.03 

5.25 

5.46 

0.03 

4.26 

4.51 

0.04 

5.16 

5.45 

0.03 

4.22 

4.32 

0.02 

5.41 

5.63 

0.03 

5.10 

5.36 

0.03 

5.29 

5.55 

0.02 

5.30 

5.50 

0.03 

5.18 

5.44 

0.03 

5.28 

5.51 

0.03 

5.12 

5.29 

0.03 

5.31 

5.52 

0.03 

5.24 

5.46 

0.04 

5.19 

5.46 

0.03 

5.32 

5.52 

0.03 

5.25 

5.36 

0.04 

5.28 

5.53 

0.07 

4.99 

5.24 
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Table  3.  Magnitudes  of  Soviet  PNE's 

Event 

#  of  Signals 

Magnitudes  [fTtza] 

(Date) 


AZG22APR66 


AZG01JUL68 


AZG22DEC71 


AZG25APR75 


AZG29JUL76 


AZG30SEP77 


AZG170CT78 


AZG180EC78 


AZG17JAN79 


AZG14JUL79 


AZG240CT79 


Ns  Nn  Nc 


3  10  0 


43  10  3 


12  0  2 


1  16  0 


41  5  7 


21  30  1 


10  0  4 


10  0  1 


S.E.M. 


0.05 


0.03 


0.05 


0.05 


0.03 


0.03 


0.06 


0.06 


0.05 


0.06 


0.06 


0RN220CT71 

31 

7 

5 

0.03 

4.88 

ORN30SEP73 

25 

9 

3 

0.03 

4.82 

ORN10JUL83A 

25 

12 

1 

0.03 

4.92 

ORN10JUL83B 

28 

10 

0 

0.03 

4.90 

ORN10JUL83C 

23 

13 

2 

0.03 

4.82 

ORN21JUL84A 

7 

4 

1 

0.06 

4.92 

ORN21JUL84B 

7 

4 

1 

0.06 

4.83 

ORN21JUL84C 

7 

4 

1 

0.06 

4.87 

PNE21MAY68 

41 

9 

1 

0.03 

4.93 

PNE29AUG74 

27 

18 

0 

0.03 

4.27 

AZG:  Azgir;  ORN:  Orenburg. 


Table  3.  Magnitudes  of  US  Explosions  outside  NTS 

Event 

#  of  Signals 

Magnitudes  [/TTzg] 

(Date) 

Ns  Nn  Nc 

S.E.M. 

Pa 

Pb 

P 

•^max 

CANNIKIN 

48  0  20 

0.02 

6.52 

6.78 

7.01 

LONGSHOT 

70  4  3 

5.09 

5.45 

5.81 

MILROW 

53  0  19 

0.02 

6.07 

6.32 

6.61 

FAULTLESS 

47  0  3 

0.03 

6.07 

6.40 

6.69 

GASBUGGY 

11  36  0 

0.03 

4.46 

4.64 

4.90 

RIOBLANCO 

15  20  0 

0.03 

4.27 

4.74 

5.02 

RULISON 

9  36  0 

0.03 

4.29 

4.41 

4.77 

SALMON 

6  33  0 

0.03 

3.85 

4.32 

4.56 

SHOAL 

16  27  0 

0.03 

4.62 

4.78 

5.04 
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Table  3.  Magnitudes  of  French,  Indian, 


S.E.M. 


0.05 


0.03 


Event 

#Of 

Signals 

(Date) 

Ns 

Nn  Nc 

BERYL 

11 

6 

0 

CORUNDON 

11 

39 

0 

EMERAUDE 

14 

23 

0 

GRENAT 

32 

30 

1 

OPALE 

3 

48 

0 

RUBIS 

45 

4 

0 

SAPHIR 

55 

3 

5 

TOURMALINE 

27 

37 

0 

TURQUOISE 

11 

51 

0 

TU19FEB77 


TU19MAR77 


TU24NOV77 


TU30NOV78 


TU25JUL79 


TU23MAR80 


TU19JUL80 


TU03DEC80 


TU25JUL82 


TU19APR83 


TU25MAY83 


RAJ18MAY74 


and  Chinese  Explosions 


Magnitudes  [/nz  g] 


Pb 


4.99 


16  25  0 


20  4  1 


31  0  0 


18  0  0 


27  12  3 


37  1  2 


31  9  0 


22  12  0 


20  1  0 


17  0  0 


7  23  0 


CH22SEP69 

27  15  0 

CH270CT75 

12  24  0 

CH170CT76 

12  33  0 

CH140CT78 

16  32  0 

CH04MAY83 

2  33  0 

CH06OCT83 

16  12  1 

CH03OCT84 

10  12  0 

CH19DEC84 

3  110 

CH05JUN87 

19  3  12 

CH29SEP88 

2  24  0 

4.60 

4.90 

5.24 

4.47 

4.65 

4.84 

4.38 

4.48 

4.78 

4.43 

4.45 

4.86 

4.39 

4.33 

4.42 

4.91 

5.16 

5.37 

4.66 

4.88 

5.14 

4.36 

4.32 

4.56 

5.72 

5.99 

6.21 

4.49 

4.56 

4.61 
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4.6  Receiver  and  Path  Effects  on  P  Waves  from  Sem4>alatlnsk  and  Novaya  Zemlya 

Along  with  the  event  values,  the  station  terms  and  the  path  terms  are  generated  by  “migim"  at 
one  single  inversion.  These  path  and  station  corrections  are  printed  in  ASCII  format,  and  can  be  con¬ 
verted  to  a  tabulated  form  easily.  Table  4  lists  the  WWSSN  station  corrections  and  path  corrections  for 
explosions  in  nine  Eurasian  nudear  test  sites.  Note  that  the  station  terms  are  appficable  to  other  source 
regions  of  the  world  as  well.  Applying  these  path  and  station  corrections  to  any  individual  explosions 
would  yield  a  reduction  in  the  fluctuationai  variation  of  station  magnitudes  with  a  factor  ranging  from  1 .2 
to  3.  Most  Novaya  Zemlya  events  have  a  typical  reduction  factor  of  2. 

Figure  1  shows  our  receiver  terms  which  are  interred  jointly  along  with  the  source-size  estimates 
and  path  terms  from  the  worldwide  explosions.  'Hte  receiver  corrections  derived  with  our  approach 
match  the  average  tectonic  structure  underneath  e^  station  very  well,  mainly  due  to  the  broad  cover¬ 
age  of  azimuths  at  each  station.  Generally  speaking,  the  station  terms  are  positive  in  shield  regions 
such  as  Australia,  Canada,  India,  and  Scandinavia,  and  they  are  negative  in  the  east  Africa  rift  valleys, 
mid-ocean  ridges  {e.g. ,  Iceland  and  Azores  Islands),  island  arcs  {e.g. ,  Indonesia,  Japan,  and  Taiwan), 
and  Himalaya  Mountain  Ranges  (Chaman  Fault,  northern  India,  Nepal,  and  Burman  Arc).  Solomon  and 
Toksoz  (1970)  and  many  other  studies  {e.g. ,  Evemden  and  Clark,  1970;  Booth  et  al. ,  1974)  found  that 
for  stations  in  U.S.,  the  attenuation  is  higher  between  the  Rockies  and  Cascades,  and  in  the 
northeastern  U.S.  This  pattern  is  also  observable  in  Figure  1  (see  also  North,  1977).  As  North  (1977) 
put  it,  it  is  gratifying  that  a  simple  parameter  such  as  mi,  can  be  utilized  to  reveal  the  tectonics.  It 
should  be  noted,  however,  that  our  empirical  station  terms  also  include  the  effect  due  to  the  crustal 
amplification  if  such  local  site  effect  is  shared  by  all  ray  paths  from  different  test  sites  to  a  particular  sta¬ 
tion.  This  could  be  the  reason  of  a  few  outliers  such  as  HNR  (Honiara,  Solomon  Islands),  PMQ  (Port 
Moresby,  East  Papua  New  Guinea),  RAB  (Rabaul,  New  Britain),  and  BAG  (Baguio  City,  Luzon,  Philip¬ 
pines)  which  do  not  show  negative  station  terms  as  would  be  expected  from  the  strong  seismicity  in  that 
region  (cf.  Figure  1).  Another  possible  reason  is  that  these  stations  have  relatively  poorer  azimuthal 
sampling  in  our  data  set,  and  hence  the  station  bias  at  these  three  stations  is  not  well  constrained.  The 
minor  discrepancy  between  the  deterministic  corrections  by  Marshall  et  at.  (1979)  and  our  empirical 
corrections  could  be  due  to  the  same  reason. 
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MEAN  STATION  AMPUFICATION  ON  mb 
252  events  used  In  ML4  Ifiversion 

1671 6+1 0055-I-2004  pati'  3,  T^2  stations  (each  recorded  10  signals  or  mors) 
ML4:  Joint  Inversion  of  source  size,  receiver  term  &  path  effect 
Assuming  each  raw  mb(ij)  s  E(i)  +  S(J)  +  F(k(l)J)  +  error 
=s>  S(J)  [receiver  term]  s  network  mean  of  [mb(iJ)-E(l)-F(k(l)J)] 

Polar  azimuthal  equidistant  projection,  78.00 , 50.00 


JunM1993  ^ 
UmrIHi 

SW  SMign:  Jlh  11/91 


1)  the  station  bias  which  needs  to  be  correcMd  (in  addition  to  the  path  eHect). 


2)  BSW  -  SW  subsite,  Balapan;  BNE  -  ME  subsite,  Balapan;  BTZ  •  transition  zona.  Balapan;  Deg  -  Degeien  Mountain;  Mzk  -  Murzhik;  NNZ  -  northern  island.  No- 


vaya  Zemlya;  Azg  -  Azgir,  Om  «  Orenburg:  PRC  >  Lop  Nor. 
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1)  the  station  bias  which  needs  to  be  corrected  (in  addition  to  the  path  effect). 


2)  BSW  -  SW  subsite,  Balapan:  8NE  -  NE  subsite,  Balapan;  BTZ  -  transition  zone,  Balapan;  Deg  -  Degelen  Mountain;  Mzk  -  Murzhilc  NNZ  -  northern  island,  No- 


vaya  Zemlya;  Azg  -  Azgir,  Om  -  Orenburg;  PRC  -  Lop  Nor. 
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1)  the  itatjon  bias  which  needs  to  be  corrected  (in  addition  to  the  path  enact). 


2)  BSW  •  SW  subsite,  Balapan;  BNE  -  ME  subsito,  Balapan;  BTZ  -  tiansltion  zone,  Balapan;  Deg  •  Degelan  Mountain;  Mzk  -  Murzhik;  NNZ  -  northern  island,  No- 


vaya  Zemiya;  Azg  -  Azgir,  Om  «  Orenburg;  PRC  -  Lop  Nor. 
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Figures  2  through  7  show  the  map  of  the  “pure  path  effect  '  (top)  and  the  combined  station 
arTH>lification  (bottom)  (defined  as  the  sum  of  the  receiver  term  and  the  path  effect)  for  explosions 
detonated  in  six  source  regions  in  Northern  Novaya  Zemiya  and  Eastern  Kazakhstan  wNch  irK:lude 
Degelen  Mountain  [Deg]  and  Mutzhik  [Mzk]  in  addition  to  the  three  subregions  defined  in  Ringdai  et  al. 
(1992):  southwestern  Balapan  [BSW],  northeastern  Balapan  [BNE],  and  the  transition  zone  [BTZ] 
between  BSW  and  BNE.  The  path  term  at  each  stsfiion  can  be  regarded  as  the  azimuthal  variation 
(towards  the  various  source  regions)  relative  to  the  averaged  station  amplification.  An  important  obser¬ 
vation  is  that  all  these  five  test  sites  exhibit  very  dfiferent  azimuthal  and  radial  amplitude  variations. 
Events  from  Degelen,  Murzhik,  and  BTZ  are  systematically  enhanced  in  the  western  U.S.  and  reduced 
in  the  eastern  U.S.,  whereas  events  from  BSW  and  BNE  are  enhanced  in  essentially  the  whole  of  U.S. 
Murzhik  events  are  reduced  in  Scandinavia,  but  Balapan  and  Degelen  events  get  strongly  enhanced 
there.  Such  highly  directiorKfependent,  distance-dependent,  and  site-dependent  patterns  of  the 
amplitude  fluctuation  could  be  a  diagnostic  for  the  path  effects  in  the  proximity  of  the  test  sites.  Back 
projections  [e.g. ,  Lynnes  and  Lay,  1990)  of  the  ntf,  residuals  onto  the  upper  mantle  and  the  lower  crust 
reveal  that  similar  mi,  residuals  come  into  alignment  in  several  regions  partitioned  by  known  geological 
features  (Jih  and  Wagner,  1991a).  Murzhik  events  recorded  in  the  western  U.S.  and  in  northeast  Asia, 
Degelen  events  in  the  western  U.S.,  and  SW  Balapan  events  at  western  European  stations  must 
pass  through  the  area  between  Chinrau  fault  and  Chingiz-Kalba  shear  zone.  All  these  paths  show 
positive  mi,  residuals.  The  area  north  of  Chinrau  fault  might  have  some  complex  features  that  result  in 
negative  mean  mi,  residuals.  Paths  from  NE  Balapan  to  North  America  and  many  continental  Euro¬ 
pean  stations  must  cross  this  area  or  even  travel  along  the  Chinrau  fault  before  entering  the  deeper 
mantle,  and  hence  the  complexity  in  the  waveforms  is  inevitable.  It  seems  that  the  mean  -Lg 
separation  of  0.07-0.17  m.u.  (e.g.,  Ringdai  and  Hokiand,  1987;  Ringdai  et  al.,  1992;  Richards  et  al., 
1990;  Section  4.6  of  this  report)  between  the  NE  and  SW  subregions  of  Balapan  could  be  due  in  part  to 
the  path  effects,  in  addition  to  the  difference  of  source  medium  postulated  previously  by  Marshall  et  al. 
(1984).  A  detailed  discussion  on  the  seismic  variabiUty  within  Balapan  test  site  is  given  in  a  later  sec¬ 
tion.  Path  effects  can  also  explain  why  the  SW  Balapan  waveforms  tend  to  be  more  complex  at  YKA 
than  those  recorded  at  WRA,  EKA,  and  GBA  arrays. 

The  initial  P  waves  from  the  three  adjacent  test  sites  have  virtually  the  same  incident  angle  at 
each  teleseismic  station,  and  anything  in  common  aooss  all  events  (such  as  the  cnjstal  amplification  as 
well  as  the  upper  mantle  attenuation  underneath  the  receiver)  would  have  been  lumped  into  the  constant 
station  term.  Thus  the  station  residuals  averaged  over  all  events  from  the  same  test  site  would  correlate 
very  little  with  the  receiver.  Instead,  they  should  reveal  more  site-dependent  information  about  the 
focusing/defocusing  pattern  underneath  E.  Kazakhstan. 

The  largest  and  most  prominent  fault  in  the  region  is  the  southeast-trending  Chingiz  right-lateral 
strike-slip  fault  that  passes  about  10  km  southwest  of  Degelen  Mountain  and  right  across  the  Murzhik 
test  area  (Rodean,  1979;  Bonham  et  al. ,  1980;  Leith,  1987b).  Soviets  reported  that  this  fault  has  a  very 
steep  dip,  which  is  consistent  with  its  linear  expression  over  large  distance  as  seen  on  Landsat  imagery 
(Bonham  et  al.,  1980).  A  distinct  fault-line  scarp  is  developed  along  much  of  the  oldest  metamorphic 
rocks.  Chingiz  Fault  extends  for  a  total  length  of  about  700  km.  Soviet  reports  postulate  that  this  fault 
extends  down  to  the  boundary  of  the  granite  layer  of  the  crust  and  possibly  into  the  upper  mantle.  For 
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Murzhik  explosions,  the  propagation  of  P„  and  Lg  waves  could  be  affected  by  this  fault  significantly, 
which  results  in  a  radiation  pattern  such  as  we  observe  in  Figure  3.  More  specifically,  the  rays  towards 
the  NW  direction  could  be  reflected  or  diffracted  to  other  quadrants,  due  to  its  post-criticai  incidence 
angles.  Such  relatively  distant  crustal  stmcture  should  have  little  impact  on  the  first  P  waves  of  Balapan 
explosions  at  teleseismic  distances,  however.  As  a  result,  amplitudes  of  Balapan  events  recorded  at 
Scandinavian  stations  are  still  largely  controlled  by  the  weak-attenuating  shield  paths  {cf.  Figures  4,  5, 
and  6). 

Marshall  et  al.  (1992)  analyze  Oegelen  and  Murzhik  events  recorded  at  4  U.K.-designed  arrays, 
and  they  find  that  EKA  and  GBA  have  distinguishable  path  effects  for  these  two  test  sites.  Amplitudes 
of  Murzhik  events  are  significantly  reduced  at  EKA,  whereas  those  of  Degelen  events  are  magnified. 
On  the  other  hand,  GBA  shows  a  strong  enhancement  for  Murzhik  signals,  but  nearly  no  effect  on 
Degelen  events.  At  YKA  or  WRA,  the  station/path  effects  are  about  the  same  for  Oegelen  and  Murzh3( 
explosions.  All  these  observations  (Figure  6  of  Marshall  at  at.,  1992)  are  in  excellent  agreement  with 
our  result  based  on  WWSSN  recordings.  The  following  is  excerpted  from  Table  4,  which  illustrates  the 
distinct  path  effects  at  EKA  and  GBA.  Note  that  the  consistent  trend  across  stations  of  wide  spatiai 
spread  as  illustrated  in  Tables  5  and  6  suggests  that  these  path  effects  are  due  to  some  very  near¬ 
source  focusing/defocusing  feature. 


Table  5.  Path  Terms  for  Stations  Close  to  EKA 

Station 

Path  Correction 

A' 

Code 

Degelen 

Murzhik 

(km) 

ESK 

0.092 

-0.525 

3.4 

VAL 

0.093 

-0.130 

602 

KON 

0.064 

-0.315 

903 

COP 

0.053 

-0.504 

985 

*)  A:  distanoe  from  EKA. 


Table  6.  Path  Terms  for  Stations  Close  to  GBA 

Station 

Path  Correction 

A 

Code 

Oegelen 

Murzhik 

(km) 

KOD 

0.021 

0.323 

373 

POO 

0.119 

0.239 

666 

NDI 

0.045 

0.155 

1669 

*)  A:  distance  from  GBA. 


The  inferred  path  terms  for  Novaya  Zemlya  explosions  have  been  compared  against  the  travel¬ 
time  residuals  to  characterize  the  propagation  paths  (Jih  and  Wagner,  1992a).  The  results  indicate  that 
paths  from  the  northern  test  site  in  Novaya  Zemlya  to  stations  in  North  America  have  systematically  fas¬ 
ter  arrivals  and  smaller  amplitudes,  suggesting  a  profound  defocusing  effect  on  the  first  arrivals;  while 
stations  in  Ireland,  Scotland,  Spain,  Bangladesh,  northern  India,  Pakistan,  Korea,  and  Kenya  report  slow 
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arrivals  and  large  amplitudes,  suggesting  a  focusing  effect.  Amplitudes  for  paths  to  Greenland.  Iceland, 
Alaska,  Turkey,  Germany,  Luzon,  Zimbabwe,  Italy,  Pueto  Rico,  Ethiopia,  and  Hawaii,  however,  seem  to 
be  controlled  by  the  anelastic  attenuation  with  slow  rays  also  associated  with  small  amplitudes,  and  fast 
rays  associated  with  large  amplitudes. 
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STATION  AMPLIFICATION  OF  mb  FOR  DEGELEN  SHOTS 

Top:  spatial  pattern  of  Isolated  path  effect 
Bottom:  receiver  term  +  path  effect 
252  events  used  In  ML4  Inversion 

1671 6+1 0055+2004  paths,  132  stations  (each  recorded  10  signals  or  more) 
ML4:  )oint  inversion  of  source  size,  receiver  term  &  path  effect 
Assuming  each  raw  mb(ij)  s  E(i)  +  SO)  +  F(k(0,J)  +  error 
Polar  azimuthal  equidistant  projection,  78.00 , 50.00 
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STATION  AMPLIFICATION  OF  mb  FOR  MURZHIK  SHOTS 

Top:  spatial  pattern  of  isolated  path  effect 
Bottom:  receiver  term  path  effect 
252  events  used  in  ML4  inversion 

16716+10055+2004  paths,  132  stations  (each  recorded  10  signais  or  more) 
ML4:  Joint  inversion  of  source  size,  receiver  term  &  path  effect 
Assuming  each  raw  mb(i,J)  =  E(0  +  S(j)  +  F(k(i)J)  +  error 
Polar  azimuthal  equidistant  projection,  78.00 , 50.00 
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STATION  AMPLIFICATION  OF  mb  FOR  BSW  SHOTS 

Top:  spatial  pattern  of  Isolated  path  effect 
Bottom:  receiver  term  +  path  effect 
252  events  used  In  ML4  Inversion 

16716+10055+2004  paths,  132  stations  (each  recorded  10  signals  or  more) 
ML4:  Joint  Inversion  of  source  size,  receiver  term  &  path  effect 
Assuming  each  raw  mb(l,J)  =  E(l)  +  S(])  +  F(k(l),])  +  error 
Polar  azimuthal  equidistant  projection,  78.00 , 50.00 
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Path+Receiver  Path  Effect 
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STATION  AMPLIFICATION  OF  mb  FOR  BTZ  SHOTS 

Top:  spatial  pattern  of  isolated  path  effect 
Bottom:  receiver  term  +  path  effect 
252  events  used  in  ML4  Inversion 

167164-10055+2004  paths,  132  stations  (each  recorded  10  signals  or  more) 
ML4:  Joint  inversion  of  source  size,  receiver  term  &  path  effect 
Assuming  each  raw  mb(l,J)  s  E(i)  +  S(|)  +  F(k(i)J)  +  error 
Polar  azimuthal  equidistant  projection,  78.00 , 50.00 
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Path+Receiver  Path  Effect 
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STATION  AMPLIFICATION  OF  mb  FOR  NNZ  SHOTS 

Top:  spatial  pattern  of  Isolated  path  effect 
Bottom:  receiver  term  4'  path  effect 
252  events  used  In  ML4  Inversion 

1 671 6-I-1 00554-2004  paths,  132  stations  (each  recorded  10  signals  or  more) 
ML4:  Joint  Inversion  of  source  size,  receiver  term  &  path  effect 
Assuming  each  raw  mb(IJ)  =  E(l)  4  SO)  4  F(k(l)J)  4  error 
Polar  azimuthal  equidistant  projection,  55.00 , 73.50 
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4.7  Plotting  Maps  of  Station  and  Path  Tenns;  gaomap 

With  the  exception  of  Figures  9  and  12.  all  other  figures  in  this  report  are  generated  by  the  same 
plotting  routine.  “Geomap"  is  a  simple  routine  to  plot  symbols,  faults,  rocks,  and  uncertairty  ellipses  on 
the  maps.  It  also  reads  (x,y,z)  pairs  from  standard  input  and  generates  PostScript  codes  to  draw  sym¬ 
bols  at  (x.y)  with  size  scaled  by  z.  The  positive  and  negative  z  data  are  drawn  in  crosses  and  octagons, 
respectively.  The  program  also  superimposes  the  plot  on  a  map  which  includes  cunres  (rivers,  faults, 
boundaries),  other  symbols,  labels,  and/or  rocks  (polygons).  Some  typical  command  calls  of  tNs  routine 
look  like 

geomap  [Ql-a  AneJl<  OaeJl-e  EtHeJi-g  GOe] H LUell-m  Mae](-p  Pf»e] 

[-S  StUe}<  input(x.y,((x,y)] 

geomap  [-fj  [-area  anaJD]  (-label  labels]  (-map  map]  (-ellpse  atipses]  (-pmj  projection] 

(-gray  blobs]  (-symb  symbols]  (-curve  curves]  <  inptti(x.y.HKy)] 

geomap  (-IJ  (-a  areaJD]  (4  labels]  (-m  map]  (-p  projection]  (-etiipse  etipses] 

(-pattern  bhbs]  (-s  symbols]  (<  curves]  (-marry  Nttel  IBei ...  HeN] 


All  the  arguments  are  optional,  and  they  are  insensitive  to  the  order.  The  first  2  sample  command 
lines  shown  above  are  for  the  case  of  one  single  dsrta  file  (/.a. ,  one  map).  The  third  is  used  for  several 
sets  of  data  to  be  plotted  on  separate  figures  (with  the  same  options)  superimposed  on  the  same  map. 
There  are  several  auxiliary  inputs  to  specify  the  options: 

Ane:  regional  ID  label  on  tire  figure 

CfUe:  curves  to  be  drawn  on  tire  figure 

BfOe:  eKpses  to  be  drawn  on  tire  figure 

GfUe:  polygortal  blabs  to  be  shaded  or  to  be  titied  witir  predefined  patterns 

Lfile:  extra  ASCII  labels  urtdemealh  tire  figure 

Mfite:  map  to  be  superimposed  on  the  figure 

Pftte:  projection  parameters  (boundary,  center  etc) 

Sfile:  extra  symbols  to  be  ptottod  on  the  figure 


-a  or  -area:  The  label  for  different  areas  on  the  map  has  a  format  very  similar  to  that  for  -a  except  in  5 
columns:  x  &  y  (at  which  the  string  starts),  ASCII  region  ID,  size  (in  inches),  and  the  direction  of  the 
label  (in  degrees).  For  example: 


77.56 

50.06 

"MurzIulC  0.09 

0.00 

78.05 

49.65 

"Degelerf 

0.09 

0.00 

78.87 

49.8 

"Balapaif 

0.09 

0.00 

77.52 

49.9 

"Chingiz  Strike-Slip  FaulT 

0.08 

■45.0 

77.86 

50.03 

"Chingiz-Kalba  Shear  Zotrd" 

0.08 

-30.00 

78.31 

50.17 

"Chinrau  FaulT 

0.08 

-25.00 

•c  or  -curve:  The  curve  (fault/river)  file  is  a  concateneftion  of  arbitrarily  many  segments  with  the  same 
format.  Each  segment  has  one  line  specifying  the  number  of  points  in  this  segment  as  well  as  the  thick¬ 
ness  of  the  curve/fault,  which  is  then  followed  by  (x.y)  pairs.  For  example, 

8  s 

77.6518  49.7423 
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77.B4S2 

49.5577 

77.9196 

49.4962 

77.9643 

49.4423 

78.0636 

49.3806 

78.1429 

49.3346 

783214 

49.2654 

784107 

49.2500 

6 

5 

77.4821 

50.2500 

77.5179 

50.2038 

77.6071 

50.1423 

77.6518 

50.1115 

77.7857 

50.0038 

77.9018 

49.9423 

•e  or  -ellipse:  The  ellipse  file  has  7  columns  defined  as  follows:  [1]  &  [2]  X.Y:  center  of  ttte  elipse.  wtth 
the  same  unit  as  the  data  on  the  map.  These  could  thus  be  in  degrees  or  km  or  whatever.  [3]  (firection 
(in  degrees):  the  direction  along  which  the  sennimaiof  axis  will  be  rotated.  [4]  &  [5]  semi-maior  and 
semi-minor  axes  (in  inches).  [6]  pen  number  (integer):  the  pen  code  (for  the  boundary  of  the  eWpse) 
runs  from  0  through  11,  the  same  as  those  hi  llbposLa.  [7]  gray  code:  an  integer  between  0  and  255 
which  determines  the  gray  level  inside  the  ellipse  (0— >  darkest,  255—>brfghtest).  The  centers  of  the 
ellipses  will  be  transformed  to  the  desired  ooordhiate  system  under  which  the  map  and  data  are  plotted 
out.  However,  the  ellipse  itself  will  not  be  transformed.  A  sample  uncertahty  ellipse  file  looks  like: 


77.700 

50.00 

*30.0 

0.40 

0.20 

2 

1 

78100 

49.75 

■45.0 

0.20 

0.10 

9 

255 

78950 

49.95 

*15.0 

0.30 

0.15 

6 
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*g>  -gray,  or  -shadow:  The  blob  file  is  a  concatenation  of  arbitrarily  many  segments.  Each  segment  has 
one  line  specifying  the  number  of  points  in  this  segment  as  well  as  the  gray  level  of  the  polygon,  which 
is  then  followed  by  (x,y)  pairs.  For  example. 


6 

0.1 

77.9018 

49.9423 

77.9643 

49.8962 

781875 

49.8192 

78.3393 

49.7885 

78.4554 

49.7654 

78.5000 

49.7500 

-pattern:  This  option  serves  the  same  purpose  as  does  -g  except  that  the  polygons  are  filled  with 
selected  patterns.  The  format  is  the  same  except  that  the  gray  level  is  replaced  by  any  integer  between 
1  and  36.  Each  of  these  codes  represents  a  predefined  pattern. 

-many  :  if  turned  on,  the  program  expects  to  read  multiple  data  buffers  so  that  a  map  win  be  drawn  for 
each  data  set  with  the  same  background  settings. 

-I  or  -label:  The  program  will  read  extra  labels  for  the  whole  plot.  The  first  line  specify  how  many  lines 
of  labels  will  be  printed  out  at  the  bottom  of  the  figure,  which  is  followed  by  ASCII  character  strings  as 
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specified.  For  example; 


6 

MEAN  STA  TION  AMPUFICA  VON  ON  mb 
253  events  used  in  MLS  inversion 

16716+10055+2004  paths,  132  stations  (each  recorded  10  signals  or  more) 
MLS:  joint  inversion  of  source  size,  receiver  term  4  path  effect 
Assuming  each  raw  mb(ij)  =  E(i)  +  S(j)  +  F(k(i),j)  +  error 
==>  S(j)  (receiver  term]  =  network  mean  of  [mb(i,i)-£(i)-F(k(i).j)J 


-m  or  'map;  The  map  file  is  a  concatenation  of  arbitrarily  many  segments.  Each  segment  has  one  line 
specifying  the  number  of  points  in  this  segment,  which  is  followed  by  (x,y)  pairs  with  format  {8f9.3).  For 
example, 


22 

65.000 

25.315 

64.324 

25.338 

64.724 

25.348 

64.592 

25.148 

64.504 

25.243 

64.263 

25.246 

64.118 

25.368 

63.864 

25.338 

63.622 

25.363 

63.476 

25.238 

63.236 

25.208 

62.971 

25.248 

62.713 

25.266 

62.519 

25.256 

62.365 

25.188 

62.244 

25.138 

62.111 

25.203 

61.968 

25.098 

61.825 

25.088 

61.725 

25.038 

61.631 

25. 136 

61.615 

25.158 

74.463 

A 

36.959 

74.450 

37.074 

74.770 

37.268 

4 

73.711 

36.906 

74.031 

36.835 

74.307 

36.897 

74.483 

36.959 

•p  or  -proj:  Currently  there  are  6  projection  methods  installed: 

1  ==>  Linear  projection, 

2  ==>  S^reographic  projection  centered  at  a  given  point 

3  ==>  Polar  azimuthal  equidistant  projection  , 

4  3s>  Far-apex  conical  projection  centered  at  a  given  point 

5  ==>  radial  plot 

6  ==>  McCartor  projection. 

The  first  line  of  the  projection  file  always  specifies  the  selected  method.  A  sample  file  for  projection 
option  1  (and/or  6)  looks  like; 

1 

77.0  79.6  49.4  50.5 

The  2nd  line  gives  the  coordinates  of  the  bottom-left  and  top-right  corners:  Xmin,  Ymin,  Xmax,  and 
Ymax.  For  other  projection,  the  2nd  line  gives  the  coordinate  of  the  projection  center  as  well  as  the 
radius  (in  degrees)  of  the  map  in  degrees.  Option  4  will  show  only  up  to  90  degrees  by  definition.  A  typi¬ 
cal  example  would  be: 

3 

76  50  100.000 


-s  or  -symb:  The  symbol  file  has  4  columns  representing  x,  y,  symbol  code,  and  size  (in  inch)  respec¬ 
tively.  The  symbol  code  runs  from  1  through  19,  the  same  as  those  in  libpost.a,  plus  a  few  extra  sym¬ 
bols:  in  particular,  100  =  filled  circle,  101  =  filled  star,  -100  =  blank  circle,  and  -101  =  blank  star.  A 
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sample  symbol  file,  “Sfile",  looks  Hce: 


79.0 

5500 

101 

0.00 

Etas 

-116.4 

37.25 

101 

500 

MTS 

179.0 

51.0 

101 

om 

AmehOka 

47.8 

450 

101 

500 

Atgt 

53.3 

51.4 

101 

500 

Onnburg 

550 

755 

101 

500 

MHZ 

5.05 

24.0 

101 

0.00 

SMtan 

-130.0 

-22.4 

101 

500 

Tuamoki 

71.7 

250 

101 

500 

Inda 

853 

41.4 

101 

aos 

LapMcr 

•f  or  -fun:  (for  projection  1  and  6  only).  H  turned  on,  the  x  and  y  axes  wM  have  different  scale. 

•z  or  -bouiMl:  The  (positive)  value  provided  after  this  flag  is  used  as  ttw  maximum  z  value  to  determtoe 
the  scale  of  symbol  size.  Normal  scale  setting  is  0.35‘/max. 

Figure  1  is  a  typical  example  of  ptotting  the  GLM  station  terms.  The  station  amplifications  and  the 
path  corrections  (tor  each  spedfic  source  regior^  are  pidted  on  top  of  the  world  map  with  the  following 
script: 


gtomp  -m  WORLD  -p  PR»  -tURt  -»  SRa  -z0.5<  Rev_Bhet>  A  Ener 

where  “Lfile”,  “Sfile",  and  “Pfile"  are  the  label,  extra  symbols,  and  the  projection  method,  respectively, 
given  In  the  discussion  above.  The  Input  data  file.  *‘Rcv_Effect”.  Ists  the  ooortftiate  and  receiver  term 
of  each  station  (cf .  Tdiie  4).  The  size  of  the  receiver  corrections  is  normalized  by  a  preset  value  of  0.5 
(cf.  the  argument  “-z  0.500”  in  the  command  line). 


-171.777  -13.909 

0.100 

AFI 

-70.4150  -23.705 

0.070 

AMT 

-71.4010  -16.462 

0.146 

ARE 

-122.235  *37.877 

0.077 

arcs 

-147.793  *64.900 

-.006 

COL 

-956020  *30.479 

0.054 

JCT 

-53.5330  *69.250 

-.121 

GDH 

-90.3000  -0.7330 

0.047 

aiE 

-105.371  *39.700 

-.236 

QOL 

*144.912  *13.SM 

-.224 

auA 

(lines  dslelad) 

- 
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4.8  Comparison  of  Various  Magnitudes 


Bocharov  et  al.  (1989)  released  the  source  information  of  96  historical  Soviet  nuclear  tests  con¬ 
ducted  in  Central  Asia  during  1965-1972.  Their  list  was  promptly  translated  arid  published  in 
EOS,  Trans.  A.G.U.  by  Vergino  (1989).  The  following  tables  are  adapted  from  those  of  Vergino's  with 
our  /T]b(PmM)  appended  as  the  column  “TG". 


Table  7.  Shagan  River  (Baiapan)  Region 

Date 

Lat 

Depth 

Yield 

Rock 

ISC 

NEIS 

Sykes 

UK 

TG 

(N) 

(E) 

(m) 

(kt) 

nh> 

/TJp 

mp 

/np 

mp 

650115 

49.9350 

178 

100-150 

Sa 

5.8 

6.3 

5.905 

5.931 

5.90 

78.9855 

316 

<20 

Sa 

5.350 

5.354 

5.30 

691130 

78.9558 

472 

125 

Co 

5.954 

6.048 

5.95 

710630 

49.9460 

78.9805 

217 

<20 

Co 

5.2 

WBM 

5.290 

5.027 

5.04 

720210 

50.0243 

78.8781 

295 

16 

Al 

Bl 

5.370 

5.370 

5.35 

721102 

49.9270 

78.8173 

521 

165 

Al 

6.1 

6.2 

6.181 

6.224 

6.16 

721210 

50.0270 

78.9956 

478 

140 

TS 

6.0 

6.0 

5.989 

5.996 

6.03 

Table  8.  Konystan 

[Murzhic)  Region 

Lat 

Long 

Depth 

Yield 

Rock 

ISC 

NEIS 

Sykes 

UK 

TG 

(N) 

(E) 

(m) 

(kt) 

nib 

mb 

mb 

nib 

mb 

651014 

49.9906 

77.6357 

048 

1.1 

mm 

661218 

49.9246 

77.7472 

427 

20-150 

Po 

5.800 

5.922 

5.88 

670916 

49.9372 

77.7281 

230 

<20 

Sa 

5.3 

5.300 

5.245 

5.21 

670922 

49.9596 

77.6911 

229 

10 

Al 

^1 

5.3 

5.200 

5.160 

5.15 

671122 

49.9419 

77.6868 

227 

<20 

A) 

4.8 

4.800 

4.410 

4.38 

681021 

78.4863 

31 

Ar 

681112 

49.7124 

78.4613 

31 

0.2x3 

Gs 

690531 

49.9503 

77.6942 

258 

<20 

A) 

5.3 

Bfl 

5.300 

5.14 

691228 

49.9373 

77.7142 

388 

40 

Al 

Bl 

wsm 

5.700 

5.791 

5.78 

700721 

49.9524 

77.6729 

225 

<20 

Sa 

Bl 

5.4 

5.400 

5.376 

5.31 

701104 

49.9892 

77.7624 

249 

<20 

Po 

Bl 

5.4 

5.400 

5.439 

5.38 

710606 

49.9754 

77.6603 

299 

16 

Al 

■-iM 

5.480 

5.526 

5.45 

iKTiTiCT 

49.9690 

77.6408 

290 

<20 

Po 

Bl 

KOI 

5.410 

5.538 

5.42 

711009 

49.9779 

237 

12 

Al 

5.3 

5.4 

5.320 

5.371 

5.25 

711021 

77.5973 

324 

23 

Sa 

5.6 

5.510 

5.580 

5.47 

720826 

49.9820 

77.7166 

285 

<20 

Al 

5.3 

1^31 

5.370 

5.363 

5.29 

720902 

49.9594 

77.6409 

185 

2 

Sa 

KOI 

5.1 

4.880 

4.788 

4.71 

Gr  Granite,  OP  CSuartz  Porphyrite,  Sa  -  Sandstone,  Al  -  Aleuroiite  (Sfltstone) 
Po  <  Porphyrite,  OS  »  Quartz  Syenite,  Gs  -  Gritstone,  Ar  =  Argillite  (Mudstone) 
Co  =  Conglomerate,  TS  =  Tuffaceous  Sandstone 


Table  9.  Degelen  Mour^nous 


611011 


640315 


640516 


640719 


641116 


650303 


650511 


650617 


650729 


650917 


651008 


651121 


651224 


660213 


660320 


660421 


660507 


660629 


660805 


660819 


660907 


661019 


661203 


670130 


670226 


670325 


670420 


670528 


670629 


670715 


670804 


671017 


671030 


671208 


680107 


680424 


680611 


Lat 


(N) 


49.77272 


49.77747 


49.81597 


49.80772 


49.80908 


49.80872 


49.82472 


49.77022 


49.82836 


49.77972 


49.81158 


49.82592 


49.81919 


49.80450 


49.76164 


49.80967 


49.74286 


49.83442 


49.73667 


49.76431 


49.82708 


49.82883 


49.74711 


49.76744 


49.74569 


49.75361 


49.74161 


49.75642 


49.81669 


49.83592 


49.76028 


49.79436 


49.81714 


49.75442 


49.84519 


49.79300 


77.99500 

116 

78.00164 

238 

78.07517 

220 

78.10197 

253 

78.09292 

168 

78.13344 

194 

78.05267 

196 

77.99428 

103 

78.06686 

152 
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Tables. 


Lat 


(N) 


49.75469 


49.82264 


49.74161 


49.81197 


49.80053 


49.74594 


49.82147 


49.75942 


49.74603 


49.81564 


49.77631 


49.78250 


49.73367 


49.79558 


49.74781 


49.73131 


49.80150 


49.80972 


49.75975 


49.74564 


49.79847 


49.76853 


49.80164 


49.74342 


49.82639 


49.76003 


49.74531 


49.73306 


49.82675 


49.73750 


49.76547 


49.81939 


49.73919 


(E) 


78.08994 


78.07447 


78.07558 


78.12194 


78.13911 


78.09203 


78.06267 


78.07578 


78.11133 


78.12961 


77.99669 


78.09831 


78.10225 


78.12389 


77.99897 


78.09861 


78.10681 


78.12839 


78.00539 


78.09917 


78.10897 


78.03392 


78.13883 


78.07850 


77.99731 


78.03714 


78.11969 


78.07569 


78.11547 


78.11006 


78.05883 


78.05822 


78.10625 


Degelen 


Depth 


(m) 
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208 
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290 
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701217 


710322 


710425 


710525 


711129 


711215 


711230 


720310 


720607 


720706 


720816 


721210 


721228 


Gr  >  Granite,  QP  »  Quartz  Porphyrita,  Sa  =  Sandstone.  Al  c  Aleurolite  (SUtstone) 
Po  °  Porphyrite,  QS  ^  Quartz  Syenite,  Gs  >  Gritstone,  Ar  =  Argillite  (Mudstone) 
Co  =  Conglomerate,  TS  =  Tuffaceous  Sandstone 
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Table  10  compares  and  Pt,  relative  to  at  several  Eurasian  nudear  test  sites.  Note  that 
there  appears  to  be  a  bias  of  0.10  m.u.  in  mtiPma^-f^biPa)  b^een  Eastern  Kazadth  and  Novaya 
Zemlya.  This  bias  couM  be  largely  due  to  the  difference  in  pP  interference  at  these  two  test  sites  (Jih 
and  Wagner,  1992ab). 


Table  10.  /n(,(Pmax)  and  m^iPi,)  vs.  /r^(P«)  (with  only) 

Test  Site 

mb{Pb)-mt,{P,) 

nib(^ma*)  ni^(Pj) 

« 

BSW 

0.271  ±0.006 

0.49110.006 

48 

BNE 

0.23510.023 

0.43110.031 

19 

BT2 

0.30210.017 

0.51310.029 

10 

Deg 

0.28710.012 

0.51310.014 

21 

Mzk 

0.29810.017 

0.52810.019 

13 

KTS 

0.27410.006 

0.49110.006 

111 

NNZ 

0.21810.010 

0.39210.010 

30 

Om 

0.16810.010 

0.42610.021 

8 

Azg 

0.41010.049 

0.68610.058 

10 

PRC 

0.16210.043 

0.406±0.063 

13 

Nuttli  (1987,  1988)  suggests  that  there  is  a  bias  of  about  0.2  m.u.  between  Degeien  and 
Balapan,  with  Degeien  explosions  having  even  larger  mi,  excitation  (relative  to  Lg ).  We  do  not  see  such 
Degelen-Balapan  bias  with  Nuttli’s  mi,{Lg)  (Taiide  11)  or  RMS  Lg  measured  at  NORSAR  (Table  12). 
The  Degeien  data  set  alone  is  too  small  for  decisive  conclusion.  However,  if  we  treat  MurzhOc  as  part  of 
Degeien,  as  did  Nuttli  (1987),  the  average  mi,{P„^-RMS  (NORSAR)  bias  between  Degeien  and 
Balapan  is  only  0.02  m.u.,  which  is  insignificant. 


Table  11.  m2.g  vs.  r7te,(L0)  (Nuttli)  at  Various  Sites 

Test  Site 

mb{Pa)-mb(Lg) . » 

inb(Pb)-mb{Lg) , » 

^b{P mmd  ~l^b{Lg)  t  # 

BSW 

-0.51310.023  28 

-0.23710.020  28 

-0.02510.019  28 

BNE 

-0.47810.045  14 

-0.22510.042  15 

-0.02910.036  15 

BTZ 

-0.47510.039  6 

-0.19110.031  6 

0.01510.026  6 

Deg 

-0.50810.124  5 

-0.182+0.112  5 

0.06310.099  5 

KTS 

-0.49910.019  53 

-0.22310.018  54 

-0.01410.016  54 

NNZ 

-0.56010.032  25 

-0.34210.036  25 

-0.16710.033  25 

4.9  Variability  within  Balapan  Test  Site 

Marshall  et  al.  (1984)  found  that  explosions  in  the  northeast  and  southwest  portioiis  of  Balapan 
test  site  produce  distinctly  different  waveforms  when  recorded  at  the  UK  seismologicai  array  stations, 
suggesting  that  Balapan  test  site  can  be  subdivided  into  two  areas  characterized  by  different  geophysi¬ 
cal  properties.  Ringdal  and  Hokland  (1987)  find  that  this  pattern  is  persistently  present  whether 
based  on  worldwide  network  or  mu  (Pcoda)  of  NORSAR  is  used.  They  inferred  the  average  mt,  -Lg 
between  SW  and  NE  subregtons  as  0.17  m.u.  In  a  follow-up  study.  Ringdal  and  Fyen  (1988)  suggest 
that  there  appears  to  be  a  transition  zone  between  the  NE  and  SW  subregbns.  Ringdal  et  al.  (1992) 
recomputed  the  SW-NE  bias  as  0.15  m.u.  with  101  Balapan  events  recorded  at  ISC  stations  and  NOR¬ 
SAR.  Although  Ringdal  et  al.  (1992)  agree  that  the  possibility  of  a  mb{Lg)  bias  contributing  to  this 
difference  between  SW  and  NE  cannot  be  entirely  ruled  out,  they  propose  an  empirical  approach  to 
correct  for  this  bias  by  assuming  this  bias  is  solely  (be  to  a  relative  mb  bias  between  these  two  areas. 

We  followed  the  zoning  of  Ringdal  et  al.  (1992)  in  partitioning  Balapan  test  site  into  three  regions: 
southwest  (SW),  transition  zone  (TZ),  and  northeast  (NE).  Figure  8  shows  the  spatial  pattern  of  m^  -Lg 
residuals  of  Semipalatinsk  explosions  based  on  Geotech’s  values  and  RMS  Lg  values  reported  at 
NORSAR.  There  is  a  significant  difference  in  the  source  medium  across  the  Chinrau  fault  separating  the 
northeastern  and  southwestern  portion  of  Balapan  test  site,  as  reported  by  Ringdal  et  al.  (1992)  and 
Marshall  et  al.  (1984)  as  noted  in  Table  12.  The  mean  bias  between  SW  and  NE  Balapan  is 

about  0.07  m.u.  Figure  8  also  indicates  that  SW  events  near  the  edge  of  the  test  site  tend  to  have 
larger  Lg  excitation  (and  hence  negative  m^-Lg  residual).  Although  this  seems  to  be  reasonable,  we 
must  be  cautious  as  this  interpretation  is  highly  dependent  on  the  accuracy  of  the  location  as  well  as  the 
geological  information. 

Note  that  the  mt,{Pmui)-Lo  bias  of  0.07  m.u.  between  SW  and  NE  (cf.  Table  12  and  Figure  9)  is 
significantly  smaller  than  that  of  previous  studies.  Regressing  the  RMS  Lg  furnished  by  Israelson 
(1992)  and  our  m2.g  on  the  yields  published  by  Brx^harov  et  al.  (1989)  (and  Vergino,  1989)  shows  that 
NE  explosions  have  positive  Lg  residuals  and  negative  residuals,  whereas  SW  explosions  show  the 
opposite  trend  (Figure  10).  A  three-dimensional  geological  model  of  the  Balapan  test  site  by  Leith  and 
Unger  (1989)  shows  a  distinct  difference  between  the  NE  and  SW  portions  of  the  test  site,  with  the 
granites  closer  to  the  surface  and  the  alluvium  thinner  in  the  southwest.  The  thk^er  alluvium  layer  in  NE 
region  could  increase  the  waveform  complexity  and  reduce  the  magnitudes  measured  with  P^ax  •  The 
first  motion  should  be  least  affected  by  this  factor,  however.  We  suggest  that  the  mb  -Lg  bias  between 
SW  and  NE  Balapan  can  be  tentatively  decomposed  into  several  parts: 

[I]  Difference  in  pP  between  SW  and  NE, 

[II]  Difference  in  m^,  coupling,  i.e. ,  mj,  (SW)  >  m^  (NE), 

[IIIJ  Difference  in  Lg  coupling,  i.e.,  Lg  (NE)  >  Lg  (SW), 

[IV]  Effects  due  to  the  station-station  correlation  stmcture. 

[V]  Effects  due  to  the  uneven  geographical  clustering  of  stations,  as  well  as  any  path  effect  which  is  not 
fully  accounted  for  through  the  network  averaging. 

Based  on  our  /hag ,  [I]  is  about  0.03  m.u.  (cf.  Table  3).  whereas  [II]  and  [III]  are  about  0.02-0.03  m.u. 
each  (Figure  3).  The  bias  of  0.07  m.u.  for  mb(Pmax)  (Table  2)  is  essentially  the  sum  of  [I]  through  [III], 
It  reduces  to  0.03  if  based  on  the  first  motion  is  used  (Table  12). 
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For  iSC  data,  we  estimate  that  [V]  is  about  0.02  m.u.  if  /n2.2  derived  by  the  oonventioruil  LSMF 
are  used.  When  m2.g  is  used,  this  term  is  eliminated,  and  hence  a  smaller  -L^  bias  is  obtained.  [11] 
and  [III]  can  be  easily  illustrated  with  regressions  on  Bocharov’s  published  yields,  as  explained  earlier 
(see  also  Figure  10).  There  are  only  a  handful  of  Balapan  events  with  published  yields  in  Bocharov 
et  al.  (1989).  However,  the  5  large  historical  events  (for  which  the  yields  were  exchanged  ckjring  JVE) 
can  also  provide  some  supplementary  due  in  support  of  our  postulated  hypotheses  [I]  through  [III].  The 
yield  estimate  based  on  tor  two  (out  of  three)  historical  events  in  SW  subregbn  (790804B  and 
791223B)  is  larger  than  that  based  on  P,.  On  the  other  hand,  the  two  events  in  NE  subregion 
(791028B  and  840526B)  have  a  smaller  yield  estimate  based  on  Pmax  as  compared  to  P, .  The  larger 
bias  of  0.15  m.u.  that  Ringdal  et  al.  (1992)  obtained  with  m,,  (ISC)  could  have  been  sightly  “enhanced" 
due  to  [IV]  and  [V].  The  mu  determination  procedure  presented  in  this  study  does  not  correct  for  [IV] 
either.  However,  the  ct^rrtribution  of  inter-station  correlation  alone  is  believed  to  be  insignificant  if 
WWSSN  is  used. 

In  Figure  1 1  we  show  the  difference  of  path  effeds  between  BSW  and  BNE  at  each  WWSSN  sta¬ 
tion,  which  is  a  measure  of  the  relative  bias  between  BSW  and  BNE  along  each  path.  Positive  symbols 
represent  the  stations  where  BSW  events  are  enhanced  relative  to  BNE  events.  If  the  raw  station  mag¬ 
nitudes  are  used  in  the  network  averaging  without  fully  accounting  for  such  path-effect  differential, 
significatt  bias  (relative  to  the  Lg  magnitude)  will  be  present.  ISC  network  is  dominated  by  western 
European  stations,  and  hence  the  effect  due  to  [V]  woukf  be  more  severe  than  that  on  WWSSN. 


Table  12.  /n2,g  vs.  RMS  Lg  (NORSAR)’  at  Various  Sites 

Site 

m,,{P,)-mb{Lg),  IP 

( ^max)  ^ 

BSW 

-0.473±0.008  42 

0.01310.009  42 

BNE 

-0.49910.028  15 

-0.25910.024  16 

-0.05610.015  16 

BTZ 

-0.22910.016  8 

-0.02510.013  8 

Deg 

-0.46910.046  5 

-0.19410.042  5 

0.02410.034  5 

Mzk 

-0.01910.032  3 

KTS 

-0.48610.009  73 

-0.00710.007  74 

NNZ 

-0.52710.019  15 

-0.12810.023  15 

1)  from  Ringdal  and  Fyen  (1991)  and  Ringdal  et  al.  (1992).  2)  f:  number  of  events. 
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AVERAGED  SW-NE  BIAS  AT  WWSSN  STATIONS 
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Poiar  azimuthai  equidistant  projection,  78.00 , 50.00 
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5.  m^- YIELD  REGRESSION  WITH  UNCERTAIN  DATA:  dwisq  (dOlsqS) 

The  standard  approach  of  yield  estimation  is  to  use  known-yield  events  to  construct  a  maonitude- 
yield  relationship  which  is  then  utilized  to  estimate  the  yield  of  other  events.  Typically  either  the  yield  or 
the  mb  is  assumed  to  be  exact  in  the  regression.  In  reality,  however,  both  the  yields  and  the  magni¬ 
tudes  are  subject  to  error.  The  regression  result  could  be  misleading  if  we  simply  assume  that  the 
yields  of  19  Semipalatinsk  explosions  published  in  Soviet  literature  are  exact.  It  has  been  speculated 
that  Soviets  might  have  rounded  8  of  the  announced  19  yields  to  the  nearest  5  kt  or  10  kt.  An 
announced  yield  of  100  kt  (e.g. ,  660320D  in  Table  9)  could  mean  something  actually  measured  between 
95  kt  and  104  kt.  It  could  also  indicate  that  perhaps  100  kt  was  the  designed  energy  release,  and  the 
actual  yield  was  somewhere  nearby.  Likewise,  the  “real  yield”  of  2  kt  (e.g. ,  720902M  in  Table  9)  could 
be  something  between  1.5  kt  and  2.4  kt.  Below  100  kt,  the  rounding  errors  could  overwhelm  the 
presumed  standard  measurement  error  —  assuming  the  announced  yields  are  not  otherwise  “fudged". 

A  more  general  regression  routine  is  given  in  this  section  to  take  the  rounding  arvf  standard  errors 
in  the  yields  into  account.  For  each  (m^ ,  yield)  pair,  we  use  a  rarxlom  number  generator  to  produce  a 
perturbed  (mb  •  yi^kl)  pair  according  to  their  uncertainty  distribution.  A  standard  least-squared  regression 
is  then  performed  for  each  data  set  of  19  perturbed  pseudo-observations.  The  procedure  is  repeated  for 
several  hundred  iterations,  and  all  the  resulting  cal^ation  cunres  are  then  used  to  infer  the  ensemble 
behavtor.  This  “doubly-weighted  least-squares  scheme"  [DWLSQ]  is  an  extension  to  the  “ordinary 
weighted  least-squares”  [OWLS]  in  which  only  errors  in  the  rr^  would  be  used  to  adjust  the  inferred 
parameters. 

The  “upper  95%  confidence  limit"  of  the  predicted  mb  at  a  given  lcg(yieid)  level  (say,  Yq)  can  be 
computed  as  follows: 

rhb(max)  -i-  t(D.O.F..0.975)Ia2(mb )  +  o2(regression)(-l  +  ^)]®®  .  [19] 

2rf(Yj  -  Y) 

where  N  =  number  of  data  points  used  in  the  regression,  D.O.F.  N-2,  o(mb )  -  the  mean  S.E.  in  the 
network  m^  used  in  the  regression,  a(regression)  «  the  a  of  residuals,  ilib(max)  -  estimate  of  the  larg¬ 
est  possible  mean  m^  at  the  given  log(yield)  level,  Y  is  the  mean  log(yield)  used  in  the  regression,  and 
t(D.O.F.,  0.975)  is  the  97.5  percentile  of  Student’s  t  distribution  at  “D.O.F."  degrees  of  freedom.  Most 
statistics  textbooks  have  a  table  of  such  values  after  Fisher  and  Yates  (1963).  Several  commonly 
quoted  t(D.O.F.,  0.975)  values  are  listed  in  Table  13. 


Table  13. 97.5  Percentile  of  t  Distribution 

D.O.F. 

5 

10 

20 

30 

40 

60 

oo 

t(D.O.F.) 

2.571 

2.228 

2.086 

2.042 

2.021 

2.000 

1.960 

The  “lower  95%-confidence  limit"  can  be  computed  in  a  similar  way: 

mb(min)  -  t(D.O.F..0.975)[a2(mb )  +  o2(regression)(^  +  ^  )] 


[20] 
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In  the  case  where  each  error  range  in  the  X  and  the  Y  is  reduced  to  zero,  then  all  the  random 
resamplings  will  simply  produce  identical  replica  of  the  original  data  set.  Consequerttly,  the  several  hun¬ 
dred  regressions  will  all  give  the  result  identical  to  a  single  call  of  the  starxJard  least-squares.  This  illus¬ 
trates  how  the  “doubly  weighted  least-squares"  [OWLSQ]  would  degenerate  to  the  standard  least- 
squares  when  the  uncertainties  in  the  X  and  Y  shrink  to  zero.  By  the  same  reasorwig,  it  is  also  an 
extension  to  the  conventional  “weighted  least-squares”  in  which  only  the  errors  in  the  Y's  would  be  used 
to  adjust  the  inferred  parameters. 

The  “doubly  weighted  least  squares  regression  routine"  is  bnplemerted  as  “doIsqS"  under  this  pro¬ 
ject.  A  typical  user  command  looks  like  the  following: 

doIsqS  [-b  or  -j]  [-a  Add]  l-z  »_of_cycles}  (  x  »_of_cydesJ  <  IN  >  A  errorjnsg 

All  the  arguments  in  the  command  line  are  optional,  and  they  are  insensitive  to  the  order.  If  neither  -b 
nor  -i  is  provided,  the  regression  reduces  to  the  corwentioruil  least  squares.  The  flag  -b  turns  on  the 
Monte-Carlo  iteration  for  both  X  and  Y.  It  implies  If  the  flag  *1  Is  given,  but  not  -b,  then  the  resam¬ 
pling  is  conducted  for  the  yields  only.  In  this  case,  the  uncertainty  in  magnitudes  will  still  be  used  in 
weighting  the  observation  matrix.  However,  the  central  values  of  magnitude  will  not  be  perturbed. 

-a:  additional  data  to  be  plotted  as  a  reference.  No  effect  on  the  regression.  The  additnnal  data  file 
has  the  same  format  as  the  input  file. 

•x:  forcing  the  range  of  X's  to  cover  so  many  “log  unit  cycles".  The  program  will  draw  a  plot  with  X 
axis  covers  the  minimal  X  plus  so  many  magnitude  unit. 

-z:  forcing  the  range  of  Y’s  to  cover  so  many  “unit  cycles".  The  program  will  draw  a  plot  with  Y 
axis  covers  the  minimal  Y  plus  so  many  magnitude  unit. 

The  following  sample  input  file  includes  the  19  Semipalatinsk  explosions  for  which  the  yields  were 
published  by  Bocharov  et  al.  (1989)  (cf .  the  2nd  column).  The  “0"  in  the  4th  column  indicates  that  the 
error  in  the  3rd  column  is  Gaussian,  which  is  assumed  to  be  10%  is  this  case.  If  this  flag  is  “1",  the 
error  in  the  2nd  column  would  represent  the  rounding  error.  The  5th  and  6th  columns  are  our  m,,  (Pmax) 
and  the  associated  error  {cf.  Table  3). 


"6S1121D" 

29.00 

2900 

0 

5.4640 

0.0240 

“660213D- 

125.00  12.500 

0 

6. 1620 

0.0230 

••6603200- 

100.00  10.000 

0 

5.9270 

0.0230 

••670922M" 

10.00 

1.000 

0 

5.1460 

0.0230 

-6809290- 

60.00 

6.000 

0 

5.7210 

0.0240 

•6907230- 

16.00 

1.600 

0 

5.2600 

0.0240 

••691130B- 

125.00  12.500 

0 

5.9520 

0.0270 

••69122BM- 

46.00 

4.600 

0 

5.7790 

0.0250 

••7104250- 

90.00 

9.000 

0 

5.9030 

0.0290 

••710606M" 

16.00 

1.600 

0 

5.4520 

0.0260 

••711009M- 

12.00 

1.200 

0 

5.2530 

0.0290 

••711021M- 

23.00 

2.300 

0 

5.4690 

0.0300 

••720210B" 

16.00 

1.600 

0 

5.3460 

0.0290 

••7203280- 

6.00 

0.600 

0 

5.0650 

0.0280 

••7208160- 

8.00 

0.800 

0 

5.0050 

0.0280 

••720902M" 

2.00 

0.200 

0 

4.7120 

0.0290 
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•72110^  t6S.00  18.SO0  0  6.1840  0.02S0 

“721210Br  140.00  14.000  0  60340  0M2S0 

‘080014B"  110.00  11.000  0  6.0480  00320 

The  regression  result  wtth  flag  4)  on  is  shown  in  Figure  12.  In  this  sampie  njn,  we  have  aiso  turned  on 

the  flag  hi  to  include  a  dummy  data  point,  which  is  considered  as  an  outfier. 

•V60S07D-  4.00  0.400  0  4.5450  0.0320 


54 


5.1  YWd  Estimates  of  Semipalatinsfc  Explosions 

It  is  fortuitous  to  have  the  source  information  released  by  Bocharov  et  at.  (1989)  (and  Vergirw, 
1989)  to  calibrate  the  Semipalatinsk  test  site.  The  small  scatter  around  the  foOowing  cali)ration  curves 
based  on  the  regression  of  our  path-corrected  tUt,  on  the  published  yields  illustrates  how  good  the  fit 
can  be  at  the  Central  Asian  test  site. 

=  0.794(±0.020)  log(W)  +  3.868(±0.030)  .  [21] 

/ni,(P6)  =0.796(10.020)  log(W) +  4.158(10.032)  .  [22] 

/nfc(^fnax)  =  0.764(10.019)  log(W)  +  4.426(10.031)  .  (23) 

Figure  12  shows  the  regression  of  /n2.a(PnMx)  on  the  the  Soviet  yields  published  by  Bocharov 
et  at.  (1989),  which  correspond  to  Equation  [23].  The  uncertainties  in  the  mt,  s  and  the  yields  are  taken 
into  account  through  800  bootstrap  resamplings.  The  darkened  bundle  is  actually  the  collection  of  all 
800  regressions,  each  produced  by  a  possbie  realization  of  19  perturbed  (rr^ ,  yield)  pairs.  The  95% 
confidence  band  (shown  as  2  curves  around  the  darkened  bundle)  is  narrower  near  the  centroid  and 
wider  towards  both  ends,  as  expected.  The  individual  95%  confidence  intervals  of  the  two  inferred 
parameters  (/.e. ,  the  slope  and  the  intercept  of  the  calibration  curve)  are  shown  with  the  dashed  line  in 
the  scatter  plot  (bottom).  Note  that  the  dashed  rectangle  is  not  the  joint  90%  confidence  intenral,  how¬ 
ever,  due  to  the  highly  correlated  nature  of  the  two  parameters.  Degelen  event  660507D  is  not  included 
in  these  regressions,  as  suggested  by  Jih  and  Wagner  (1991, 1992b). 

We  have  utilized  these  calibration  cun/es  to  estimate  the  yield  of  aU  114  Semipalatinsk  explosions 
in  our  data  set,  and  the  results  are  summarized  in  Table  14.  For  cratering  events  (such  as  6501 15B) 
the  yield  estimate  based  on  the  first  imtion  (/.e. ,  P« )  should  be  used,  since  no  depth  correction  such  as 
that  used  in  Marshail  et  at.  (1979)  has  been  applied  to  nribiPb)  or  mk,(Pmax}  in  Table  3.  For  this  particu¬ 
lar  event,  Myasnikov  et  at.  (1970)  gave  a  "scaled  apparent  radius"  and  scaled  depth  of  of  51  and  50 
m/kt*’^,  respectively.  Ck)mbining  this  information  vwth  the  crater  radius  and  the  emplacement  depth 
released  at  the  IAEA  symposium,  Ringdal  et  at.  (1992)  inferred  the  yield  of  this  explosion  as  120  kt, 
which  is  almost  identical  to  our  estimate  of  119  kt  based  on  (Table  14).  This  example  illustrates  that 
P,  from  hard-rock  test  sites  in  a  stable  region  could  be  a  very  favorable  phase  for  the  source  size 
determination. 

Much  of  the  source  information  about  the  Soviet  JVE  explosion  (88091 4B)  has  not  been  released. 
The  "New  York  Times"  (Gordan,  1988)  states  that  the  American  and  Soviet  on-site  rrreasurements  are 
said  to  give  yields  of  1 15  kt  and  122  kt,  respectively.  If  we  substitute  the  /n2.g(Pa )  of  JVE  into  Eolation 
[6],  the  mean  yield  estimate  would  be  112  kt.  Sykes  and  Ekstrom  (1989)  gave  an  estimate  of  113  kt 
based  on  the  arithmetic  average  of  mb  and  Ms .  Priestley  et  at.  (1990)  analyzed  the  Lg  amplitudes  at 
4  seismographs  near  the  Semipalatinsk  test  range:  KSU  (Karasu),  KKL  (Karkaralinsk),  BAY  (Bayanual), 
and  TLG  (Talgar),  and  they  obtained  a  mb(Lg)  of  5.96810.02.  Murphy  et  at.  (1991)  gave  a  network- 
averaged  mb  of  6.012  with  a  standard  deviation  of  0.190  across  the  network.  They  also  derived  a 
RMS  Lg  of  5.969  using  8  stations  in  U.S.S.R.,  Norway,  and  Manchuria.  This  value  is  identical  to  that  of 
Ringdal  et  al.  (1992)  based  on  NORSAR  recordings.  It  is  worth  noting  that  all  these  seismic 
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magnitudes  give  very  consistent  yield  estimates  in  the  range  100-150  kt,  as  specified  in  the  biateral 
agreemeit  signed  by  U.S.  and  Soviet  governments  before  JVE  (Richards.  1990;  Stump,  1991). 

Thera  are  15  events  In  common  in  israelson’s  (1992)  RMS  Lg  data  set  and  our  /rtz  s  data  set  for 
which  the  Soviet-published  yields  are  available.  The  correlation  between  the  RMS  Lg  and  niz  o  residu¬ 
als  (relative  to  the  expected  magnitude  at  the  associated  yield  value)  is  very  weak  and  hence  the  combi¬ 
nation  of  these  two  methods  for  a  better  yield  estimate  is  justifiable. 

It  is  Meresting  to  note  that  three  out  of  the  five  “historical  events”  (for  which  the  yields  were 
exchanged  in  1988)  have  a  yield  of  153  kt,  based  on  our  (P, )  (Table  14).  The  remaining  two  histori¬ 
cal  events  and  the  JVE  all  have  a  yield  around  115  kt,  based  on  the  first  motion.  The  yields  would  have 
a  larger  variation  for  each  of  these  two  groups,  if  the  based  on  the  more  conventional  largest  cycle 
(i-O-  •  ^max)  was  used  instead. 
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800-bootstrap  Intercepts  mb(WWSSN,ML4,Pmax) 


mb:Yield  Relationship  at  KTS 


Yieid(KT)  (Bocharov  et  ai.,  1989) 


OWLS  (uncertain  X  &  Y):  S=0.76(0.019),  1=4.43(0.031).  19.  data  uaed. 

95%  en-or  in  mb  at  1.1 0,50,1 00,1 50KT:  0.20, 0.11, 0J)9. 0.11, 0.12, 

95%  factor  in  yield  at  1.1 0,50,1 00,1 50KT:  3 JO,  1.94, 1.73, 1  Jl,  2.10 

OWLS  (precise  X  assumed):  S=0.78(0.032).  1=4.41  (0i)51) 

Standard  LS:  8=0.77(0.030),  1=4.42(0.048) 

10%  S.E.  in  yield  assumed 

Scatter  Piot  of  Inferred  Parameters 


800-bootstrap  Slopes 

95%  confidence  Interval  of  slope:  0.76-t-/-0.041 
95%  confidence  interval  of  intercept:  4.43W-0.066 

[97.5%  quantile  of  t(1 7.  D.o.F.),  2.110,  used]  Jurmu* 

Usar  Hh 

Figure  12  SW  dnign:  tih  04^1 
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Table  14.  Yield  Estimates  of 


Epicemer 


Lon  Lat 


Semipalatinsk  Explosions 


Yield  E^imate 


Pb 


6501 15B 


6602130 


6603200 


6605070 


6610190 


6702260 


67091 6M 


671122M 


68061 9B 


6809290 


690531 M 


6907230 


6909110 


691130B 


691228M 


700721 M 


701104M 


7103220 


7104250 


710606M 


710619M 


710630B 


711009M 


711 021 M 


7112300 


72021 OB 


7203280 


7208160 


720826M 


720902M 


721102B 


721 21 OB 


7212100 


730723B 


731214B 


Yield 


Announced 


100-150 


29 


125 


100 


Event 

Date 

Site 

7504278 

8TZ 

760704B 

8SW 

761 1238 

8NE 

7612078 

8SW 

770329D 

Deg 

7705298 

8SW 

7706298 

8NE 

770730D 

Deg 

7709(^8 

8NE 

7803260 

Deg 

780422D 

Deg 

8SW 

7807058 

8SW 

780728D 

Deg 

7808298 

8NE 

7809158 

8SW 

7811048 

8NE 

7811298 

8SW 

7906238 

8SW 

7907078 

8NE 

7908048 

8SW 

7908188 

8TZ 

7910288 

8NE 

7912028 

8SW 

7912238 

8SW 

800522D 

Deg 

8SW 

8009148 

8SW 

8010128 

8NE 

8012148 

8TZ 

8NE 

8104228 

BSW 

8109138 

BTZ 

8110188 

BSW 

8111298 

BSW 

8112278 

BSW 

8204258 

BSW 

8207048 

BSW 

Table  14.  Yield  Estimates  of 


Epicenter 


Lon 


78. 


78.898 


78.947 


78.840 


78.140 


78.772 


78.849 


78.160 


78.914 


78.070 


78.170 


78.802 


78.867 


78.140 


78.968 


78.862 


78.949 


78.796 


78.845 


78.992 


78.887 


78.919 


78.994 


78.786 


78.753 


78.082 


78.798 


78.811 


79.022 


78.917 


78.978 


78.807 


78.894 


78.846 


78.847 


78.780 


78.887 


78.810 


Lat 


49.939 


49.903 


50.018 


49.944 


49.790 


49.946 


50.044 


49.770 


50.059 


49.730 


49.720 


49.913 


49.903 


49.756 


50.008 


49.928 


50.046 


49.956 


49.915 


50.039 


49.903 


49.948 


49.976 


49.910 


49.933 


49.784 


49.938 


49.931 


49.968 


49.909 


50.068 


49.899 


49.914 


49.928 


49.902 


49.933 


49.918 


49.961 


Semipalatinsk  Explosions 


Yield  Estimate 


Pb 


Yield 


Announced 


HE:  historical  events  discussed  at  U.S.-U.S.S.R.  negotiation  during  'aZ-Sa. 


Table  14.  Yield  Estimates  of  Semipalatinsk  Exptosions 


Event 


Date 


820831 B 


821205B 


821226B 


83061 2B 


831006B 


831026B 


831120B 


840219B 


840307B 


840329B 


840425B 


840526B 


840714B 


840915B 


841027B 


841202B 


841216B 


841228B 


85021 OB 


850425B 


850615B 


850630B 


850720B 


87031 2B 


87041 7B 


870620B 


870802B 


871115B 


871213B 


880213B 


880403B 


880614B 


881112B 


881217B 


890708B 


Site 

Lon 

BSW 

78.762 

BSW 

78.845 

BNE 

78.995 

BTZ 

78.897 

BSW 

78.757 

BSW 

78.824 

BNE 

79.018 

BSW 

78.744 

BNE 

78.954 

BTZ 

78.919 

BSW 

78.851 

BNE 

79.004 

BSW 

78.877 

BNE 

78.911 

BSW 

78.817 

BNE 

79.007 

BSW 

78.816 

BSW 

78.692 

BSW 

78.779 

BSW 

78.881 

BSW 

78.839 

BSW 

78.668 

BSW 

78.786 

BSW 

78.826 

BSW 

78.779 

BSW 

78.670 

BSW 

78.746 

BSW 

78.875 

BSW 

78.756 

BSW 

78.793 

BSW 

78.868 

BTZ 

78.906 

BSW 

78.749 

BNE 

78.958 

BSW 

78.823 

BNE 

78.968 

BSW 

78.923 

BSW 

78.779 

Epicenter 


Lat 


49.914 


49.928 


Yield  Estimate 


49.923 


49.925 


49.912 


50. 


49. 


50.054 


49.923 


49.936 


49.968 


49.908 


49.992 


49.906 


50.010 


49.947 


49.881 


49.898 


49.925 


49.907 


49.864 


49.948 


49.936 


49.919 


49.883 


49.937 


49.881 


49.962 


49.933 


49.907 


49.950 


50.024 


8 


50.047 


49.881 


49.868 


HE;  historical  events  discussed  at  U.S.-U.S.S.R.  negotiation  during  '67-’88. 
JVE:  Joint  Verification  Experiment 


Pb 

P  max 

8 

8 

126 

135 

38 

40 

129 

129 

71 

88 

124 

122 

19 

18 

47 

44 

21 

29 

78 

81 

81 

84 

153 

131 

153 

151 

Yield 


Announced 


125 

117 

129 

139 

150 

157 

3 

3 

112 

133 

11 

11 

57 

64 

20 

23 

5.2  AsMssment  of  m,,  Bias 

We  regressed  m2.»  values  of  21  high-coupling  NTS  tests  again^  the  announced  yields  (DOE. 
1990),  and  the  resulting  calibration  curves  are  listed  betow: 

m6(P,)  (NTS)  =  0.758(±0.015)log(W) -1-3.636(10.033)  .  124) 

"»t,(Pi,)  (NTS)  =  0.825(10.015)  log(W)-h  3.699(10.033)  .  (25J 

mb(Pnmi  (NTS)  =  0.81 1(10.015)  log(W)  +  3.977(10.033)  .  [26] 

The  KTS-NTS  mt,  bias  can  then  be  computed  in  a  straightforward  manner  by  comparing  Equa¬ 
tions  [21]-[23]  against  [24]-[26].  Between  100  and  150  kt,  the  KTS-NTS  bias  is  estimated  as  0.35 

m.u.  using  our  mb(Pmax)  (cf.  Table  15).  For  comparison,  we  have  also  mcluded  in  Table  15  the  bias 

estimate  inferred  from  Murphy's  (1990, 1981)  calibration  cunres  based  on  the  network-averaged  spectra 
[NAS].  The  inversion  algorithm  Murphy  ef  al.  (1989)  adopted  in  their  NAS  scheme  is  the  conventional 
LSMF  (Equation  [13]).  The  NAS  method,  by  its  frequency-domain  nature,  excludes  clipped  or  noisy  sig¬ 
nals  in  the  magnitude  computation,  which  is  quite  different  from  the  time-domain  approaches  (such  as 
ours)  in  which  the  maximum-litelihood  method  can  be  applied  to  count  for  the  censoring  effects. 
Murphy's  (1990, 1981)  formulae  are 

mj,  (NAS)  (KTS)  =  0.75  log(W) -f  4.45  ,  and  [27] 

mj,  (NAS)  (NTS)  =  0.81  log(W)  -t-  3.92  .  [28] 

Despite  the  methodological  difference  between  the  two  techniques  it  is  very  interesting  to  note  that 
Murphy's  KTS  formula  (Equation  [27])  is  almost  identicai  to  Equation  [23].  Also,  Murphy’s  NTS  calS>ra- 
tion  curve  (Equation  [28])  has  a  slope  identical  to  that  of  Equation  [26].  There  exists  a  bias  of  0.05-0.06 
m.u.  between  [26]  and  [28],  which  causes  a  discrepancy  of  0.05  m.u.  in  our  KTS-NTS  bias  estimate  and 
that  of  Murphy's  at  150-kt  level. 


max)  (^^S) 

0.47 

m^{P^){KTS) 

0.50 

mt(P.)(KTS) 

0.37 

NAS(Murphy) 

Table  15.  Expected  mt,  Bias  Relative  to  NTS 


m22 

W22 

50KT 

100KT 

150KT 

10KT 

50KT 

100KT 

150KT 

0.44 

0.42 

0.41 

0.40 

0.37 

0.35 

0.35 

0.47 

0.46 

0.45 

0.43 

0.41 

0.40 

0.40 

0.37 

0.38 

0.38 

0.27 

0.29 

0.30 

0.31 

0.47 

0.43 

0.41 

0.40 

6.  TIME-DOMAIN  DETERMINATION  OF  L,  PATH  CORRECTION:  guessQ 

There  are  several  different  approaches  readily  available  to  determine  the  path  0(/n; 

[A]  Apply  the  coda-O  method  of  Herrmann  (1980),  as  did  Nuttli  (1988). 

[B]  Synthesize  the  path  Oohi  along  the  great-circle  path  between  the  source  and  the  receiver  using 
the  2-dimensional  map  of  that  regioa 

[C]  Apply  GLM  [General  Linear  Model]  or  LSMF  [Least  square  Matrix  Factorization]  inversion  to  infer 
the  path  corrections  along  with  the  source  terms  (Jih,  1992;  Israelson,  1992). 

Approach  [C]  would  perform  very  well  when  some  extra  reliable  information  about  the  events  {e.g. , 
the  average  of  m^  or  mb(Lg)  values)  is  available  to  constrain  the  joint  inversion  (Jih,  1992).  Here  we 
provide  another  approach  which  is  very  similar  to  [C]  exce|M  that  the  stations  are  calibrated  individually 
with  those  events  for  which  Nuttli  (1988)  already  determined  the  mt,{Lg)  values. 

In  processing  the  Lg  data  set  assembled  for  another  contract  at  Geotech,  the  “sustained  max¬ 
imum  motion"  of  Lg  phase  are  measured  in  a  manner  identic^  to  that  which  Nuttli  (1986ab,  1987, 1988) 
proposed.  That  is,  the  amplitude  equaled  or  exceeded  by  the  three  largest  amplitude  waves,  of  the 
vertical-component  Lg  waves  with  period  around  1  second  were  picked.  The  station  amplitude  reading  is 
first  corrected  for  the  effects  of  geometrical  spreading  and  dispersion  with  the  formula  appropriate  for  the 
Airy  phase;  the  residual  (relative  to  Nuttli's  mt,{Lg))  is  then  regarded  as  completely  due  to  the  anelastic 
attenuation  along  the  path: 

where  A  is  the  epicentral  distance  in  km,  A(A)  is  the  observed  Lg  amplitude  measured  in  the  time 
domain  in  pm  [microns]  at  the  epicentral  distance  of  A  km.  The  corresponding  Q(f)  can  then  be  deter¬ 
mined  in  a  straightforward  manner: 

Qffl  ■  •  [301 

where  U  is  the  group  velocity.  Once  a  suite  of  Q(f)  values  is  available  for  a  station  of  interest,  a  linear 
regression  is  then  conducted  to  find  the  maximum-likelihood  estimate  of  the  quality  factor,  Oq,  as  well  as 
the  frequency-dependency,  t\,  via  the  model: 

Q(f)  =  00*^,  or 

log[Q(f)]  =  log[Ool  +  Ti*f,  [31] 

Two  simple  FORTRAN  routines  are  cc.'ibined  to  implement  this  time-domain  calibration  procedure: 
“guessQ"  and  “domie2".  The  code  “guessO”  reads  5  groups  of  parameters: 

•o  the  source  information  which  includes  the  event  magnitude  (e.g.,  mb(Lg)),  epicenter  (latitude,  longi¬ 
tude),  and  the  event  identification  (an  ASCII  string); 

-a  the  amplitude  in  nanometers; 

-p  the  period  in  seconds; 
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•V  the  group  velocity  in  km/sec; 

■8  the  station  code. 

These  entities  can  be  interchanged  arbitrarily  In  the  command  line.  The  source  parameters  under 
the  flag  ”-o"  must  be  given  in  the  specified  order,  however.  A  sample  script  for  the  station  1ST  is  given 
below.  The  5  Novaya  Zemlya  events  recorded  at  1ST  were  rather  large  in  event  size,  and  all  of  them 
were  detonated  before  TTBT  came  into  effect. 

guessQ  -o  6.4S  73.400  54.900  66300  -m  49.5  -p  0.93  -v  3.60  -» 1ST 
guessQ  -o  6. 75  73.310  55. 140  70267  271.9  -p  1.68  -v  3.30  -*  1ST 

guessO  -o  6.68  73.380  55. 100  71270  -a  1 15.3  -p  1.48  -v3.60  -s  1ST 
guessQ  -o  6.42  73.330  55.080  72241  -a  199.6  -p  1.89  -v3.50  -s  1ST 
guessQ  -o  6.43  73.350  55.070  75294  -a  44.1  -p  0.90 -v  3.50  -s  1ST 


The  code  "guessQ"  computes  the  path  y  based  on  Equation  (29]  with  the  input  information.  The 
resulting  y  and  the  input  parameters  are  printed  out  in  a  format  as  foUows: 


Y>  0.001588  Q(f)~  590.77  />  1.08 
y.  0.00t327  Q(f)m  381.68  U  0.53 
y.  0.001505  Q(fh  391.68  f~  0.68 
y-  0.001209  Q(f)=  392.75  0.53 

y-  0.001607  0(1)=  620.60 1=1.11 


■c  6.450  73400  54.900  66300  -a  49.5  -p  0.93  -v  3.60  -s  1ST 
■o  6.750  73.310  55.140  70287  -a  271.9  -p  1.88  -v  3.30  -s  1ST 
-o  6.680  73380  55.100  71270  -a  115.3  -p  1.48  -v  3.60  -s  1ST 
-o  6.420  73.330  55.080  72241  -a  199.6  -p  1.89  -v  3.50  -s  1ST 
-o  6.430  73.350  55.070  75294  -a  44.1 -p  0.90  -v  3.50  -s  1ST 


The  columns  of  frequency,  f,  and  the  Q  are  then  extracted  from  the  output  of  “guessQ"  and  con¬ 
verted  to  the  following  form  after  taking  the  logarithm.  The  code  "domle2"  reads  this  input  file  through 
the  direct  input.  It  anticipates  to  read-in  a  free-formatted  file  which  consists  of  4  columns.  Each  line 
gives  a  data  point  to  be  used  in  the  maximum-likelihood  regression.  ■  The  second  column  is 
log(frequency),  which  can  be  measured  with  very  high  precision.  The  third  and  fourth  columns  are  the 
lower  and  upper  bounds  of  log(Q(frequency)),  respectively.  The  first  column  (in  quotes)  is  the  quality  flag 
of  Y.  Four  choices  of  this  data  quality  flag  are  permissible:  and  ”%". 


+0.0315171 

■0.2741579 

-0.1702617 

-0.2764618 

+0.0457575 


2.7714190 

2.5817049 

2.5929310 

2.5941107 

2.7928126 


2.7714190 

2.5817049 

2.5929310 

2.5941107 

2.7928126 


Running  "domle2"  to  regress  log(Q(f))  on  log(f)  assuming  a  linear  model  Y  «  A  +  B  X  ,  we  get  B  = 
0.64410.099,  A  =  2.74910.019,  a(Y)=0.032,  p=  0.9660.  Therefore,  with  the  present  procedure,  the  Lg 
path  correction  appropriate  for  Novaya  Zemlya-IST  path  is  q  =  B  -  0.644  and  Oo  =  10^  =  561  (cf. 
Table  16).  Since  all  Y  values  of  this  example  are  uncensored,  the  maximum-likelihood  estimate  of  the 
slope  and  intercept  would  be  identical  to  those  based  on  the  standard  least  squares.  A  non-trivial  exam¬ 
ple  which  involves  censored  data  is  given  in  Section  7. 
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6.1  Lg  Path  Corrections  for  Novaya  Zamlya,  Samipalatinsk,  and  NTS 

The  resulting  path  corrections  for  Novaya  Zemiya  test  site  are  listed  in  Table  16,  along  with  those 
corrections  of  Nuttli’s  (1988).  The  match  is  fairly  good.  This  simply  suggests  that  Geotech's  Lg  ampli¬ 
tude  measurements  furnished  by  Rivers  et  at.  (1993)  could  be  very  consistent  with  Nuttli's.  It  Is  krterest- 
ing  to  note  that  1ST  (Istanbul,  Turkey)  and  TP.I  (Trieste,  Italy)  did  record  Lg  phases  from  large  historical 
Novaya  Zemiya  events.  Along  with  the  7  WWSSN  stations  for  which  Nuttli  (1988)  already  published  the 
Qq  values,  now  we  have  a  total  of  12  paths  calibrated  for  Lg  waves  from  Novaya  Zemiya.  Stations 
KON  (Konsberg,  Norway)  and  KBS  (Kingsbay,  Svalbard)  are  not  well  constrained  due  to  the  limited  data 
size,  and  hence  Nuttli’s  (1986b)  Qf/r\  would  have  to  be  retained. 


Table  16.  0(/Tt  for  Novaya  Zemiya  Lg 

Station 

This  Study* 

Nuttli  (1988) 

Code 

Oo 

T1 

Oo 

i\ 

COP 

668 

0.41 

633 

0.4 

KBS 

315 

0.5 

KEV 

249 

0.74 

252 

0.6 

KON 

496 

0.5 

NUR 

433 

0.42 

420 

0.5 

STU 

550 

0.55 

531 

0.5 

UME 

397 

0.82 

391 

0.5 

ESK 

463 

0.63 

DAG 

270 

0.69 

1ST 

561 

0.64 

NOR 

223 

0.43 

TRI 

417 

0.24 

*  Based  on  ampliruds  measuremenis  furnished  by  Rivers  el  el.  (1993). 
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*  Based  on  amplituda  measuremanu  furnished  by  Rivera  at  al.  (1903). 
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*  Based  on  amplitude  measurements  furnished  by  Rivers  ef  af.  (1993). 
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7.  UNEAR  REGRESSION  WITH  CENSORED  OBSERVATIONS:  dOfnle2 

Consider  the  situation  where  the  independent  variable  can  be  precisely  measured,  and  that  we 
want  to  regress  the  dependent  variable  as  a  linear  function  of  the  independent  variable.  This  is  an 
extension  to  the  single-event  network  magnitude  determination  discussed  in  Section  1.3. 

Suppose  we  have  a  linear  rrKXfel  of  the  common  form: 

Y  =  a  +  px  +  V  ,  [321 

where  X  is  the  independent  variable  which  has  a  precision  relatively  much  better  than  that  of  Y,  the 
dependent  variable,  a  and  p  are  the  intercept  arxf  slope,  respectively,  to  be  determined,  arKf  v  is  an 
error  term,  v  is  assumed  to  be  a  Gaussian  random  variable  with  mean  zero  and  standard  deviation  o. 
Furthermore,  assume  that  there  are  four  types  of  data  available: 

[0]  the  observed  measurement,  Y,  is  known  as  yo, 

[1]  Y  is  only  known  to  be  less  than  certain  level, 

[2]  Y  is  only  known  to  be  larger  than  certain  level,  and 

[3]  Y  is  only  krrawn  to  lie  between  two  bounds. 

Type  3  data  are  not  uncommon.  The  majority  of  Soviet  yields  recently  published  by  Bocharov 
et  al.  (1989)  and  Vergino  (1989)  actually  fall  in  this  category  (cf.  Table  9). 

Elegant  maximum-likelihood  theory  can  be  derived  for  this  model.  Suppose  there  are  no,  ni,  na, 
and  no  measurements  for  each  type,  respectively.  The  conditional  IH'^'lihood  function  of  the  censored 
observations  (  yo,  ti,  ta,  to)  given  the  intercept  a,  slope  p,  and  a  is 

no  n, 

L  (  yo.  t,.  ta,  to  I  a,  P,  a )  =riP(  Yj  =  yoj  I  o,  p,  o  )  *  nP(  Yj  <  tij  I  a.  p,  o )  *  [331 

j=i  j=i 


”2  ”3 

nP(  Yj  >  taj  I  a,  p,  o )  •  nP(  taj  <  Yj  <  tbj  I  a,  P,  a ) 


i  =  i 


j  =  i 


and  the  log-likelihood  function  is 


In  L  (  yo.  ti,  ta,  to  I  a,  p,  a )  =  -  -^  ln(27ia^)  - 


■^_L(yoi-a-PXoj)^  + 


[341 


Xin  4»(zij)  +  Sin  <l>(-zaj)  +  £ln  [«I>(Zbj)-^(Zaj)l  . 

i=i  j=i  j=i 

where  Z|  =  [yj  -  a  -  pt|l/a:  y©,  t,.  ta,  and  to  are  the  vectors  of  the  four  data  types. 

Solving  s  o  implies  immediately  that  the  d.  the  optimal  estimate  of  a,  must  satisfy  the  follow- 
oo 

ing  necessary  condition: 
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^yoj-a-pxoi)* 

_ H _ _ 


Solving 


0  implies  that  the  sum  Of  the  ‘defined  residuals*’  should  be  zero.  Solving 


OVIVIliy  30t  ^  HII^IOO  UlCU  HIW  9UIII  Wl  UfW  lOffHtvu  wv  —  w 

implies  that  the  vector  of  refined  residuals  should  be  orthogonal  to  the  vectors  of  means.  It  follows  that 
the  optimal  estimate  of  a  and  p  can  be  obtained  by  the  “standard  least  squares*  inversion  wHh  the  cen¬ 
sored  data  all  replaced  by  their  conditional  expectations,  /.e. .  the  “r^ned  observations”.  Thus  o  can  be 
solved  iteratively  with  [6]  along  with  a  and  p  using  the  EM  algorithm.  In  the  non-censored  case,  this 
"dbm/e2"  code  gives  results  identical  to  those  derived  by  the  standard  least  squares. 


Example 


Lynnes  and  Baumstark  (1991)  measured  the  Pg/Lg  ratio  of  51  NTS  explosions  at  several  fre¬ 
quency  bands.  At  0.5-1 .0  Hz,  there  were  21  events  for  which  the  spectral  ratio  was  clipped,  i.e. ,  Pg  /Lg 
was  greater  than  a  certain  level,  but  it  could  not  be  precisely  detemiined.  H  only  the  30  uncensored 
spectral  ratios  were  used  in  the  least-squares  regression,  the  best  linear  fit  would  have  a  slope  of  p  - 
0.350  and  an  intercept  of  a  «  -0.111  (cf.  the  dashed  line  in  Figure  13).  This  resuit  is  misleading 
because  it  implies  that  there  is  very  little  depth  dependence  of  Pg  /Lg  ratio.  When  the  21  dips  are 
induded,  however,  the  slope  increased  signlTicanUy  (p  «  1.175,  a  >  -0.237).  This  is  a  typical  example 
illustrating  how  the  maximum-likelihood  approach  can  lead  to  a  more  reasonable  model  by  including  the 
censored  information.  In  this  case,  the  maximum-lkelihood  result  does  reveal  the  decrease  of  relative 
Lg  excitation  with  the  shot  depth  (cf.  the  black  line  ki  Figure  13). 


The  program  reads  4  columns  of  data  via  the  standard  input.  The  first  column  is  the 

data  type  in  quotes  “<”,  “>”,  or  "%").  The  second  column  is  the  independent  variable  (which  is 
assumed  to  be  precisely  known).  The  third  and  fourth  columns  give  the  measured  thresholds  of  the 
dependent  variable.  For  data  of  type  3,  the  upper  bound  arxf  the  lower  bound  are  different,  arxl  hence  2 
columns  are  required.  The  sample  input  file  shown  below  is  taken  from  Lynnes  and  Baumstark  (1991) 
(pages  19-23)  where  the  independent  and  dependeit  variables  are  the  shot  depth  and  the  iog(P0  /Lg  ) 
(measured  at  0.5-1 .0  Hz),  respectively,  of  51  NTS  explosions.  The  program  "domle2"  ignores  the 
source  information  appended  in  each  line. 

0.320  0.122  0.122  WUS  OOSIOma  3935 

"c"  0.320  0.110  0.110  WUS  81149aa  3945 

0.340 -0.3»3 -0.393  WUS  81 191  ae  3952 

0.472  0.273  0.273  WUS  81274^  3960 

"«■'  0.445  0.080  0.080  WUS  81315ac  3964 

0.494  -0.031  -0.031  WUS  81337aa  3969 

0.335-0.237-0.237  WUS  81350ad  3971 

"='■  0.400-0.126-0.126  WUS  82210ac  3975 

0.640  0.338  0.338  WUS  8221 7ad  3977 

0.229-0.066-0.066  WUS  82245aa  3979 

"="  0.406  0.095  0.095  WUS  82266ac  3981 

0.451  -0.379  -0.379  WUS  82266ad  3982 

"="  0.304  -0.075-0.075  WUS  83042aa  4036 
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0.343 -a078 -aOTS 
“m’  0.384  0.190  0.190 
0.330-0.088-0.088 
0.336-0.179-0179 
0.S33 -0104 -0104 
0.388 -0007 -0007 
0.800  0147  0147 
0.381  0.210  0210 
“m-  0.366  0.324  0324 
“m-  0.640  0317  0317 
0.S1S -0179 -0179 
“m"  0.640  0.337  0337 
0.000  02SO  02S9 
0.293  0048  0048 
0.381  -0012  -0012 
“m"  0.332 -0.124 -0124 
“m-  0.400  0113  0113 
0573  0582  0.562 
‘S'*  0.637  0.266  0JX6 
‘S"  0204  0044  0.044 
‘S"  0.294 -0078 -0.078 
‘S'*  0.213  0149  0.149 
“>"  0.518  0.309  0.309 
“>••  0.640  0.504  0.504 
“>"  0.564  0.329  0.329 
‘S“  0.366  0.210  0.210 
“>■’  0.320  0.105  0.105 
‘S-  0.625  0.301  0.301 
“>"  0.335  0.367  0.367 
••>"  0.483  0.461  0.461 
•S"  0.415  -0.091  -0.091 
‘S  '  0.579  0.462  0.462 
‘S  '  0.500  0.636  0.636 
‘S"  0.610  0.805  0.805 
‘S  ’  0.400  0.264  0.264 
‘S  '  0.600  0.562  0.562 
‘S 0.500  0.374  0.374 
‘S  '  0.500  0.646  0.646 


¥8VS03(Mab4097 
WUS83146ae4044 
¥ntS89180me4048 
WUS  83215m  4049 
¥¥US8326Sm4054 
WUS  84091m  4071 
WUS  841S2tb  4072 
WUS  84172m  4075 
WUS  a^49m  4079 
WUS  8435005  4084 
WUS85082m4087 
WUS8SO92m4O80 
WUS8S006m4091 
WUS  8516305  4005 
WUS  85177m  4097 
WUS8S229m4103 
WUS  882050  4131 
WUS00352m3038 
WUS81157ae3046 
WUS  81197m  3053 
¥nJS  8123005  3956 
WUS  8126705  3969 
WUS  81316m  3966 
Wl^  8202805  3973 
WUS82272m3986 
WUS  82316m  4018 
¥nJS  83223m  4050 
WUS  83244m  4051 
WUS  8421505  4077 
WUS  84257m  4081 
WUS  85289m  4107 
WUS  85330m  4108 
WUS85362m4111 
WUS  86081m  4114 
WUS  86100m  4115 
WUS  86112m  4117 
WUS  86141m  4118 
WUS  861  S6ad  4121 
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APPENDIX:  PREREQUISITE  MATHEMATICS  FOR  MAXIMUM-UKELIHOOD  ESTIMATOR 

Proposition  1.  Let  X  be  a  Gaussian  random  variable  with  the  probabiiity  density  function  [p.d.f.]  g  and 
the  cumulative  distribution  function  [c.d.f]  G.  respectively-  Its  mean  and  variance  are  denoted  by  p  and 
o^.  respectively.  Then 

a 

J  xg(x)dx  =  pG(a)  -  o^gCa). 


Proof. 


J  "9(x)dx  =  ^^x  exp(-  ■^)dx 


20^ 


20^  ■  V2jci. 

=  pG(a)-a2g(a). 


20^ 


In  particular,  when  a  =  oo,  this  integral  gives  the  mean  of  X,  namely,  p. 

Proposition  2.  Let  X  be  a  Gaussian  random  variable  with  mean  p  and  variance  o^,  then  E  [  X  I  X  <  a  ] 
=  p  -  a2g(a)/G(a) . 

Proof.  Let  Y  be  the  random  variable  X  I  x  <  a.  then 

P(Y<b)  =  P(X<b(X<a)  =  P(  X  <  X  <  a ) 

which  is  1  if  b  >  a,  and  G(b)/G{a)  if  b  5  a.  Therefore,  the  p.d.f.  of  Y  is: 

h(x)  =  0  if  X  >  a  ,  h(x)  =  g(x)/G{a)  if  x  <  a, 
and  the  expectation  of  Y  is 

~  a 

E(Y)  =  J  xh(x)dx  =  J  xg(x)/G(a)dx 

=  Proposition  1) 

=  p  -  a2g(a)/G(a) . 

Similarly,  it  can  be  shown  that  E[XlX>a]  =  p-i-  o^(a)/G(-a) .  Note  that  the  conditional  expec¬ 
tation  E  [  X  I  X  >  a  ]  is  the  “best”  guess  of  X  under  the  constraint  that  one  knows  only  that  X  >  a. 

In  computing  E  [  X  I  X  >  a  ],  it  is  generally  more  convenient  to  transform  the  random  variable  X 
into  the  standard  random  variable,  Z  ~  N(0,  1),  for  which  the  p.d.f.  and  c.d.f.  are  typically  available  as 
system-furnished  functions  or  part  of  some  utility  libraries  in  the  public  domain.  Let  O  and  ^  be  the  c.d.f. 

and  p.d.f.  of  Z,  respectively,  then  G(a)  =  and  g(a)  =  a  Therefore, 
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E  [  X  I  X  <  a  1  =  U- <r4(-^m-^l . 
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